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Abstract 


The  research  presented  in  this  study  is  primarily  concerned 
with  the  application  of  numerical  modelling  techniques  to 
the  stability  analysis  of  coal  mine  structures.  The  work  has 
been  limited  to  room  and  pillar  coal  mining  of  flat  and 
gently  dipping  (below  20°)  seams. 

Conventional  analytic  techniques  are  described, 
together  with  numerical  and  physical  modelling  techniques 
that  are  presently  available.  The  boundary  element  computer 
program  (DDSEAMS)  which  has  been  used  in  this  analysis  is 
described  in  detail,  and  the  results  of  limited  testing  of 
the  program  are  included.  A  stress  analysis  design  concept 
is  presented,  which  includes  the  development  of  a  stress 
criterion  for  stability  at  a  point. 

Comparisons  are  made  between  the  numerical  model,  a  3-D 
physical  model,  a  base  friction  model  and  an  underground 
phenomenological  study  at  Grande  Cache,  Alberta. 

An  illustration  of  the  stress  analysis  design  concept 
is  included,  with  reference  to  a  proposed  mine  at  Coal 
Valley,  Alberta.  This  includes  a  detailed  geotechnical 
testing  program,  and  the  analysis  of  a  number  of  alternative 
mine  geometries. 

The  major  conclusions  are: 

*  Improved  stress  analysis,  compared  to  conventional 

methods,  are  obtained  by  using  the  numerical  model 
presented . 

*  Definition  of  a  stress  criterion  is  useful  for  the 
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prediction  of  over-stressed  regions  that  are  potentially 
unstable . 

*  Where  strata  stress  conditions  or  pillar  safety  factors 

indicate  that  failure  will  take  place,  the  use  of 
elastic  analysis  techniques  is  not  adequate  for  the 
definition  of  the  overall  stability. 

*  The  3-D  physical  modelling  qualitatively  simulates  pillar 

failures  observed  in  the  underground  phenomenological 
study  at  Grande  Cache 

*  Initial  analysis  indicates  that  extraction  ratios  greater 

than  50%  may  be  possible  by  the  "punch"  mining  method 
proposed  at  Coal  Valley. 
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1.  INTRODUCTION 


The  object  of  this  research  was  to  demonstrate  the 
application  of  a  particular  numerical  analysis  technique  for 
assessing  the  stability  of  mining  structures  in  room  and 
pillar  coal  mining.  A  boundary  element  computer  program 
(DDSEAMS)  has  been  used,  which  is  simple  and  inexpensive  to 
operate,  but  which  does  have  limitations  in  its 
applicability.  These  limitations  are  not  restricted  by 
current  Knowledge,  but  by  the  time  required  to  develop 
suitable  program  packages  for  commercial  use. 

Room  and  pillar  coal  mining  is  a  method  whereby 
sections  are  developed  on  advance,  usually  consisting  of 
parallel  drives  (rooms),  and  connected  at  frequent  intervals 
to  delineate  rectangular  coal  pillars.  In  partial  extraction 
the  pillars  are  left  to  support  the  overburden.  When 
complete  extraction  is  performed,  these  pillars  are 
extracted  on  retreat,  and  the  roof  strata  are  allowed  to 
cave.  The  method  is  usually  applicable  to  flat  and  gently 
dipping  seams  (dips  less  than  20°). 

A  major  factor  in  the  choice  of  mining  method  and  mine 
layout,  is  the  stability  of  the  openings  and  pillars  that 
are  formed  during  extraction  of  the  coal.  This  study  is 
primarily  concerned  with  the  stability  of  pillars  in  room 
and  pillar  coal  mining.  However,  the  stability  of  openings 
cannot  be  ignored,  as  they  interact  with  pillars  in  mine 
structures.  Portions  of  the  analysis  will  only  be  applicable 
to  room  and  pillar  coal  mining  of  flat  and  gently  dipping 
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seams,  but  the  methods  adopted  are  applicable  to  many  other 
mining  methods  which  require  openings  and  pillars. 

A  stress  analysis  approach,  that  is  suitable  for 
solution  by  numerical  methods,  has  been  adopted.  Numerical 
problems  are  ideal  for  solution  by  digital  computers,  which 
are  rapidly  becoming  faster,  more  powerful  and  inexpensive. 
It  is  probable  that  the  future  trend  will  be  for  every  mine 
to  have  a  small  computer  facility,  or  a  communications  link 
to  a  larger  computer.  Complex  mine  stability  analyses  are 
expensive,  therefore  they  are  currently  beyond  the  finance 
and  expertise  of  small  mines. They  are  usually  restricted  to 
a  small  number  of  alternative  geometries  and  conditions.  In 
the  future  it  is  probable  that  stress  analysis  techniques 
may  be  used  as  routine  mine  design  tools.  They  enable  many 
alternative  layouts  to  be  analysed  quickly  and 
inexpensively,  with  a  minimum  of  input  data  preparation. 

Any  engineering  analysis  requires  that  certain 
simplifications  be  made  in  the  modelling  of  the  real 
situation,  in  order  to  obtain  a  solution.  Thus  the 
simulations  required  are  seldom  perfect,  and  this  is 
especially  true  for  rock  materials.  The  analysis  is  further 
complicated  by  the  limited  knowledge  of  the  true  rock  mass 
behaviour,  when  it  is  subjected  to  complex  loading 
conditions,  common  in  mining  situations.  Therefore, 
assumptions  have  to  be  made  in  defining  the  rock  mass 
behaviour,  which  are  often  further  simplified  for  analysis 
purposes.  This  produces  an  element  of  uncertainty  in  the 
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solution  that  can  only  be  reduced  by  comparison  with  the 
real  mining  conditions.  Real  behaviour  is  difficult  and 
expensive  to  determine  comprehensively.  Therefore  the 
Knowledge  of  the  pre-mining  stress  field,  and  the  stresses 
and  displacements  induced  by  mining,  are  seldom  complete. 

An  attempt  has  been  made  to  compare  the  stress  analysis 
technique  adopted  in  this  study  with  physical  models  and  an 
underground  phenomenological  study. 

A  mine  design  simulation  study  of  a  mine  property  at 
Coal  Valley,  Alberta,  was  undertaken.  A  limited  number  of 
possible  mine  geometries  were  analysed.  The  rock  mass 
parameters  were  estimated  from  geotechnical  testing  of  drill 
core  from  the  proposed  site. 


2_;_  CONVENTIONAL  DESIGN  ANALYSES 
Conventional  methods  of  analysing  the  stability  of  mine 
pillars  and  openings  will  be  defined  as  those  methods  which 
do  not  require  the  use  of  digital  computers  in  their  routine 
solution.  Numerical  analysis  techniques  which  require  the 
use  of  computers  for  their  solution  will  be  described  in  the 
next  chapter.  The  analysis  methods  have  been  divided  into 
the  methods  that  concern  pillar  stability  and  those  that 
concern  opening  stability.  This  division  is  made  for 
simplification  purposes,  but  in  reality  pillars  and  openings 
interact  in  a  complex  manner.  Only  techniques  that  are 
applicable  to  flat  and  gently  dipping  seams  (dip  less  than 
20°)  will  be  described. 

2 . 1  Pillar  Stabi 1 i tv  Analyses 

In  the  analysis  of  pillar  stability,  the  load  that  will 
be  imposed  on  the  pillar  is  estimated,  and  then  this  is 
compared  with  the  strength  of  the  pillar.  From  this 
comparison,  safety  factors  for  the  stability  of  the  pillars 
can  be  derived. 

2.1.1  Pillar  Load  Estimation 

Three  methods  for  estimating  the  loads  on  pillars  are 
presented.  These  methods  are  briefly  described: 
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Tr ibutarv  Area  Theory 

This  is  the  simplest  analytic  technique,  utilising  the 
assumption  that  each  pillar  is  loaded  by  the  vertical  column 
of  overburden  that  is  tributary  to  that  pillar.  The  load  is 
calculated  as  an  average  pillar  stress,  therefore  the  stress 
variation  within  the  pillar  is  not  accounted  for.  An  example 
calculation  for  the  square  pillar  layout  in  Figure  2.1  is 
shown  below; 


apa  =  °v ^ (Wo+wP) /Wp) 2 

o  =  Average  pillar  stress 
pa 

cv  =  Pre-mining  vertical  stress 
Wo  =  Width  of  opening 
Wp  =  Width  of  pillar 


This  equation  is  a  simple  force  balance  for  use  in  flat 
seams  where  the  vertical  normal  stress  is  a  principal  stress 
direction.  It  should  not  be  used  for  workings  that  do  not 
extend  over  sufficient  lateral  extent,  relative  to  depth 
(such  that  arching  effects  to  the  abutments  are  present),  or 
for  pillars  in  close  proximity  to  large  abutment  pillars. 

Beam  Def lect ion  Theor ies 

Beam  theories  should  be  used  when  a  roof  arching  effect 
is  present.  They  also  enable  the  non-uniform  stress 
distribution  across  the  pillar  to  be  calculated. 
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The  immediate  roof  strata  are  assumed  to  consist  of  one 
or  more  elastic  beams  (two-dimensional  problem),  which  span 
the  area  of  interest.  In  three-dimensional  analysis  the 
beams  are  replaced  by  plates,  but  the  theory  is  similar.  The 
beams  are  loaded  by  the  overburden,  which  is  assumed  to 
provide  a  uniform  load.  The  load  is  re-distributed  onto  the 
pillars,  which  are  assumed  to  be  elastic.  The  loading  of  the 
beam  is  transferred  to  the  pillars,  producing  support 
reactions.  Beam  theory  is  used  to  obtain  stresses,  by 
assuming  that  the  pillars  respond  in  a  similar  manner  to 
elastic  foundations.  Two-dimensional  analysis  is  considered 
to  be  adequate  when  the  length  of  the  beam  in  the  third 
dimension  is  greater  than  twice  the  span. 

Equations  have  been  developed  to  satisfy  different  end 
constraints,  beam  interactions  and  beam  thickness.  For  thick 
beam  analysis  ( Sheorey  &  Singh,  1974)  the  beam  thickness  is 
assumed  to  be  greater  than  one  fifth  of  the  span,  otherwise 
thin  beam  analysis  ( Stephansson ,  1971)  is  employed. 

Beam  theory  provides  a  better  understanding  of  pillar 
stress  distributions  than  the  tributary  area  theory,  and 
should  be  used  when  the  pillar  configurations  are  not 
constant.  There  are  many  assumptions  adopted  to  define  the 
roof  strata  beams  and  their  end  constraints,  which  require  a 
high  degree  of  engineering  input  and  experienced  judgement, 
to  obtain  meaningful  results.  Beam  theory  is  applicable 
where  roof  strata  beds  are  homogeneous,  and  can  withstand 
tensile  stresses.  Sedimentary  coal  measures  strata  are 


usually  weak,  bedded  and  non -homogeneous  (as  a  result  of 
jointing),  which  reduces  the  applicability  of  beam  theory. 
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E last ic  Analysis 

There  are  a  number  of  analytical  solutions  in  the 
theory  of  elasticity  which  can  be  related  to  underground 
openings.  They  are  usually  restricted  to  two-dimensional 
plane  stress  or  plane  strain  analysis.  The  material  is 
usually  assumed  to  be  an  infinite,  homogeneous,  -isotropic, 
linear  elastic  continuum.  Simple  opening  geometries,  which 
are  usually  aligned  with  the  principal  stress  field,  have 
been  solved.  Multiple  openings  can  be  simulated  by  combining 
results  from  single  openings.  The  overall  stress 
distribution  is  produced  by  summation  of  the  distributions 
for  each  individual  opening,  utilising  the  principle  of 
superposition.  This  method  can  be  used  to  produce  3-D 
solutions  by  combining  2-D  solutions  in  transverse  planes. 

Solutions  for  simple  problems  will  not  be  presented  as 
they  can  be  found  in  many  texts  on  rock  mechanics  and 
elastic  theory  (Timoshenko  &  Goodier,  1970;  Hoek  &  Brown, 
1980  and  Coates,  1967). 

2.1.2  Pillar  Strength  Cr i ter ia 

Prediction  of  pillar  strengths  is  difficult,  primarily 
because  of  the  non-uniform  nature  of  the  pillar  loading, 
resulting  in  complex  failure  mechanisms.  High  stress 
concentrations  on  the  sides  of  non-yielding  pillars  are 
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generally  encountered.  Local  failures  are  initiated  in  these 
regions  as  a  result  of  low  confining  stresses.  As  the 
perimeter  fails,  the  load  is  transferred  towards  the  central 
core,  increasing  stresses  within  this  region.  The  central 
core  is  more  confined  than  the  perimeter,  and  is  therefore 
able  to  withstand  higher  stresses  prior  to  failure.  The 
pillar  strength  is  usually  defined  as  the  maximum  stress 
that  it  can  support,  and  does  not  account  for  local  failures 
(Wilson  &  Ashwin,  1972). 

The  pillar  strength  is  a  function  of  the  material 
properties,  its  size,  shape,  end  constraints  and  loading 
conditions.  The  estimation  of  pillar  strength  is  usually 
undertaken  by  loading  small  pillars  to  failure,  and  by 
applying  size  and  shape  correction  factors  to  these  results. 

Pillar  Size  Factors 

Compressive  strength  results  are  usually  obtained  by 
testing  drill  core  samples  or  small  coal  pillars.  Core 
samples  are  usually  taken  from  the  stronger  coal  seam 
material,  because  of  the  problems  associated  with  the 
preparation  and  handling  of  the  weaker  portions.  These 
samples  have  fewer  weakness  planes  and  discontinuities  than 
the  real  pillar,  and  therefore  tend  to  have  a  higher 
strength . 

Much  work  has  been  done  in  attempting  to  define  the 
relationship  between  size  and  strength.  These  results  have 
been  summarized  by  Hustrulid  (1976),  and  a  formula  has  been 
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derived  which  seems  to  fit  the  available  data  quite  well. 
For  model  pillars  less  than  3ft.  (0.91m)  in  height: 

a  =  k//iT 
c 

ac  =  Cube  compressive  strength 

k  =  Constant  for  each  coal 

(compressive  strength  of  a  1  inch  cube) 

H  =  Pillar  height  (inches) 

For  pillars  with  dimensions  greater  than  3ft.  (0.91m),  the 
cube  compressive  strength  is  essentially  constant  and  is 
defined  by: 

°c  =  k/6 


Therefore  the  pillar  cube  strength  can  be  defined  as  one 
sixth  of  the  compressive  strength  of  a  one  inch  cube,  as 
mine  pillars  always  have  dimensions  greater  than 
3ft . (0.91m) . 

Pillar  Shape  Factors 

The  compressive  strength  of  rectangular  pillars  have 
been  investigated,  and  found  to  be  a  function  of  their 
shape.  Consequently,  the  relationships  derived  have  been 
defined  with  respect  to  the  width  and  height  of  the  pillar. 
The  majority  of  results  indicate  that  there  is  a  linear 
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relationship  between  the  pillar  compressive  strength  and  its 
width  to  height  ratio.  Hustrulid  (1975)  has  analysed  these 
results  and  found  that  one  equation  can  be  used  which  gives 
a  reasonable  fit  to  the  majority  of  available  data: 

o  =  o  (0.778  +  0.222  W/H) 

P  c 

a  =  Pillar  compressive  strength 
P 

a  =  Cube  compressive  strength 
o 

H  =  Height  of  pillar 

W  =  Width  of  pillar 

These  results  indicate  that  mine  pillar  compressive  strength 
for  a  particular  coal  is  linearly  dependent  upon  the  shape 
factor  alone,  as  the  cube  compressive  strength  of  mine 
pillars  may  often  be  assumed  to  be  constant.  Most  of  the 
studies  have  been  undertaken  on  small  samples,  which  may  not 
necessarily  reflect  the  relationship  for  the  real  mine 
pill  ars . 


2 . 2  Open i nq  Span  Stabilities 

In  room  and  pillar  mining  the  pillars  are  delineated  by 
the  openings,  consequently  the  sides  of  the  openings  are 
also  the  sides  of  the  pillars.  Therefore  this  description 
will  only  address  the  stability  of  the  roof  and  floor 
strata,  as  the  pillar  stability  has  already  been  discussed. 
Only  rectangular  openings  will  be  considered,  as  they  are 


, 


common  to  room  and  pillar  coal  mining.  If  the  span  is  found 
to  be  unstable,  then  support  must  be  installed  to  maintain 
stability. 

2.2.1  Roof  Stabi 1 i tv 

Roof  strata  can  be  analysed  by  using  beam  theory,  as 
described  previously  for  pillar  load  calculations.  The 
analysis  can  be  simplified,  as  failure  is  expected  to  be 
initiated  at  the  centre  or  abutments  of  each  span.  Therefor 
only  the  stresses  and  displacements  at  these  points  are 
required.  Expressions  for  the  maximum  normal  and  shear 
stresses  are  given  by  Duvall  (1976)  for  a  gravity  loaded 
beam  clamped  at  both  ends.  For  beams  with  a  span  to 
thickness  ratio  greater  than  five,  the  maximum  tensile 
stress  is  more  than  three  times  the  maximum  shear  stress. 
The  tensile  strength  of  rock  is  usually  much  less  than  the 
compressive  and  shear  strengths,  and  therefore  a  design 
formula  is  derived  based  on  the  tensile  strength: 

T  =  y.F.L2/2.at 

T  =  Thickness  of  beam 
y  =  Unit  weight  of  rock 

F  =  Safety  factor 

L  =  Span  length 

a  =  Tensile  strength 

The  safety  factor  in  tension  should  be  between  4  and  8 
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(Duvall,  1976).  From  this  equation  a  required  roof  beam 
thickness  can  be  calculated,  and  if  this  beam  is  not 
present,  then  a  composite  beam  should  be  formed  by  bolting 
rock  layers  together. 

Where  high  horizontal  stresses  are  present,  design  on 
the  basis  of  tensile  strength  may  not  be  applicable.  Tensile 
stresses  will  be  reduced,  and  may  even  be  eliminated.  In 
these  cases  the  design  should  be  based  on  the  shear  or 
compressive  strength  of  the  strata.  Strata  failures 
resulting  from  high  horizontal  stresses  have  been  observed 
by  Parker  (1973).  These  include  shear  failures  close  to  the 
abutments  of  the  opening,  and  compressive  failures  at  the 
centre  of  spans. 

Where  analysis  of  structural  geology  indicates  that 
jointing  may  cause  problems,  it  may  be  necessary  to  design 
supports  to  stabilise  any  wedges  formed  in  the  roof.  These 
wedges  can  be  defined  by  plotting  stereonet  projections  of 
joint  systems  and  bedding  planes. 

2.2.2  Floor  Stabi 1 i tv 

Floor  stability  problems  occur  when  weak,  thinly  bedded 
strata  are  enountered.  The  competence  of  floor  beds  is 
usually  reduced  when  water  is  present.  Floor  failures 
usually  indicate  that  high  horizontal  stresses  are  causing 
buckling  failures  of  thin  unconfined  beds.  These  failures 
are  enhanced  when  pillars  load  the  floor  strata  beyond  its 
bearing  capacity.  The  combined  failures  result  in  floor 
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migration  into  the  entry  (floor  heave).  This  is  usually  a 
time-dependent  phenomenon,  but  in  highly  stressed  regions 
closure  of  excavations  can  take  place  in  a  matter  of  days. 

Floor  heave  is  difficult  to  predict,  as  it  results  from 
both  floor  beam  failures  into  the  opening  and  floor  failures 
below  pillars.  The  interactions  between  pillar  and  floor  are 
not  well  understood,  as  complex  behavioural  mechanisms  are 
present.  As  a  result  analysis  is  usually  based  on  past 
experience  and  engineering  judgement.  Floor  heave  may  be 
controlled  by  reinforcement  of  floor  strata  (e.g.  rock 
bolting),  or  by  reducing  the  pillar  stresses  close  to  the 
opening . 


3_j.  MODELLING  TECHNIQUES 

The  simulation  of  reality  by  various  modelling  techniques  is 
discussed,  together  with  their  applicability  to  various  mine 
situations.  The  techniques  that  will  be  adopted  in  this 
study  are  described  in  more  detail. 


3 . 1  Numer ica 1  Mode 1 1 i ng  Techniques 

Numerical  modelling  techniques  have  developed  rapidly 
after  the  advent  of  the  digital  computer,  which  enables  many 
repetitive  computations  to  be  performed  on  large  quantities 
of  data.  This  allows  the  simplification  of  complex  problems 
into  a  series  of  discrete  units  which  have  Known  solutions. 
These  discrete  units  are  analytically  independent,  but  are 
connected  to  form  the  complex  problem,  and  a  solution  is 
obtained  by  satisfying  defined  boundary  conditions. 

In  analysing  physical  structures  in  rock  the  discrete 
units  may  be  boundary  elements,  finite  elements  or  blocks. 
Each  element  or  block  is  defined  as  having  a  prescribed 
relationship  between  stress  and  strain,  or  load  and 
deformation.  The  combination  of  these  units  simulates  the 
overall  mine  structure.  The  approximation  introduced  by  this 
discretisation  is  usually  adequate  for  engineering 
appl i cations . 


3.1.1  Finite  E 1 emen t  Method 

With  the  finite  element  method  the  rock  structure  is 
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divided  into  a  series  of  elements  which  are  connected  at 
nodal  points  to  form  a  framework.  In  two-dimensional 
analysis  these  elements  are  polygonal  (commonly  triangles  or 
quadrilaterals),  and  are  assumed  to  have  unit  thickness  in 
the  third  dimension.  Forces  or  displacements  applied  to  the 
element  nodes  produce  stresses  and  strains  within  the 
element,  which  are  translated  into  displacements  at  the 
nodes.  The  elements  are  connected  together  at  the  nodes  to 
form  a  continuous  structure.  The  forces  and  displacements 
must  be  compatible  at  each  node.  The  relationship  between 
forces  and  displacements  is  defined  by  a  series  of  equations 
which  depend  upon  the  material  parameters,  the  element 
geometry,  and  the  degree  to  which  the  strain  is  permitted  to 
change  across  the  element.  A  system  of  equations  is  produced 
which  is  solved  to  satisfy  the  boundary  conditions. 

The  elements  should  extend  beyond  the  area  of  interest, 
such  that  boundary  conditions  can  be  specified  that  do  not 
interact  with  the  mine  structure.  Planes  of  symmetry  are 
often  utilised  to  define  boundaries  and  reduce  the  number  of 
elements  that  are  required  to  define  the  problem  area.  The 
method  produces  continuous  variation  of  stress  and  strain 
within  elements,  but  not  between  elements.  In  order  to  model 
the  real  situation  with  sufficient  accuracy,  a  higher 
element  density  is  required  in  regions  of  rapid  stress 
variation,  such  as  excavation  boundaries. 

The  advantage  of  using  a  discrete  structure  is  that 
different  parameters  can  be  assigned  to  each  material  type, 


. 
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such  that  the  physical  structure  is  more  completely  defined. 
More  complex  rock  mass  behaviour  can  be  modelled  by  using 
special  elements  (e.g.  for  joints)  or  complex  material 
behaviour  (non-linear  elastic,  e 1  as to-p 1 ast i c ,  visco-plastic 
or  no-tension  analysis)  ( Z i enki ewi cz ,  1977). 

3.1.2  Boundary  E 1 emen t  Method 

The  boundary  element  method  involves  the  discretisation 
of  the  boundary  of  the  excavations  alone.  Integral  equation 
methods  are  used  to  produce  a  fully  continuous  solution  for 
the  stress  distribution  over  the  whole  area  (closed  form 
solution).  This  results  in  a  much  smaller  system  of 
equations  than  the  finite  element  method  and  a  consequent 
reduction  in  computer  time  and  data  preparation.  The  major 
disadvantage  of  the  boundary  element  method  is  that  it 
usually  requires  the  rock  mass  to  be  homogeneous,  and  have 
linear  elastic  properties.  It  is  possible  to  define  zones 
with  different  physical  properties  by  developing  special 
interface  elements  (Banerjee.  1976),  but  the  technique 
requires  further  development. 

3.1.3  Block  Models 

Block  models  can  be  used  to  simulate  blocky  or  jointed 
strata  as  a  discrete  system  of  rigid  blocks.  These  blocks 
interact  with  one  another  at  edges  and  corners,  where  force 
and  displacement  conditions  are  satisfied.  The  system  is 
solved  by  analysing  the  forces,  displacements,  velocities 
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and  accelerations  of  blocks  for  incremental  time  steps  until 
equilibrium  is  attained. 

Cundall  (1971)  has  developed  a  computer  model  to  solve 
a  class  of  low  stress,  strong  rock  problems;  where  block 
internal  stresses  are  not  important,  and  large  displacements 
are  to  be  expected.  In  reality,  blocks  delineated  by  joint 
systems  probably  only  contact  at  a  few  points  and  thus  the 
analysis  of  these  contact  forces  have  practical  validity. 


3 . 2  Numer ica 1  Program  ( DDSEAMS )  Adopted 

A  boundary  element  program  (DDSEAMS)  was  used  in  this 
study  (Crouch,  1976),  which  is  described  in  Appendix  1.  The 
numerical  technique  involves  using  displacement 
discontinuity  boundary  elements  and  solving  the  system  of 
influence  coefficients  for  the  defined  boundary  conditions. 
The  program  was  tested  to  ensure  that  it  was  operating 
correctly,  before  being  applied  to  any  mine  modelling. 

3.2.1  Ana  1 vt ica 1  So Tut  ion  Compar isons 

In  order  to  check  the  program  for  any  computational 
errors,  and  to  give  an  indication  of  the  accuracy  of  the 
analysis  technique,  a  simple  comparison  was  carried  out 
between  the  DDSEAMS  program  output  and  the  closed  form 
solution  for  a  thin  crack  in  a  triaxial  stress  field. 

Closed  form  solutions  for  flat  cracks  (slits),  in  an 
infinite,  homogeneous,  linear  elastic,  isotropic  continuum, 
exist  for  various  crack  geometries  and  stress  conditions. 
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The  solution  for  a  flat  elliptical  crack  was  chosen  for 
comparison;  The  long  axis  being  of  infinite  length  in  the  z 
direction  (Figure  3.1),  in  a  three-dimensional  stress  field. 
This  solution  simulates  the  plane  strain  condition,  which  is 
calculated  using  DDSEAMS  with  the  material  parameters  chosen 
to  simulate  an  isotropic  continuum.  The  analytical  solution 
takes  the  following  form  (Bray,  J.W.  unpublished  notes): 


(a)  Along  the  x  axis: 
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(c)  Along  roof  and  floor: 
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(d)  At  edge  of  crack  (sides): 
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cr  =  Normal  stress  in  x  direction 

X 

0y  =  Normal  stress  in  y  direction 

a  =  Normal  stress  in  z  direction 

=  Shear  stress  on  x  plane  in  y  di 

T^z  =  Shear  stress  on  y  plane  in  z  di 

tzx  =  Shear  stress  on  z  plane  in  x  di 

V  =  Poisson's  Ratio 


?x/  P  /  Pz  =  Field  normal  stresses 
P P-._  /  P„v  =  Field  shear  stresses 

•xy  y  z  z  x 

w  =  Width  of  slit 
(H/w  is  assumed  to  tend  to  zero) 


rection 

rection 

rection 


A  simple  program  has  been  written  to  calculate  some  of  the 
solution  results.  The  listing  of  program  (CRACK)  is 
presented  in  Appendix  2,  together  with  an  example  of  the 
output . 

A  run  of  DDSEAMS  was  made,  with  similar  conditions  to 
the  analytical  solution,  using  a  total  of  eighty  elements 
(forty  for  the  crack).  Isotropic  material  parameters  were 
used  in  the  solution,  and  the  comparison  of  results  for  the 
vertical  stress  component  (  oy  )  is  presented  graphically  in 
F i gure  3.2. 

The  most  important  stress  component  is  the  vertical 
stress,  as  this  indicates  the  pillar  and  abutment  loading. 
As  can  be  seen  from  Figure  3.2;  deviation  of  the  results 
from  program  DDSEAMS  and  the  analytical  solution,  does  not 
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exceed  10%.  The  greatest  deviation  occurs  close  to  the 
opening,  where  the  boundary  element  approximations  will  have 
greatest  effect.  Considering  the  engineering  applications  of 
program  DDSEAMS  in  mine  design,  this  deviation  is 
acceptable.  This  program  has  greatest  benefit  in  providing 
comparisons  between  alternative  mine  layouts. 

3.2.2  Program  Opt imi sat  ion 

A  number  of  similar  runs  of  program  DDSEAMS  were  made, 
simulating  identical  single  span  conditions,  but  with 
different  numbers  of  elements  defining  the  spans.  The 
results  of  these  runs  are  summarised  graphically  in 
Appendix  4. 

In  general  the  smaller  number  of  elements  that  are  used 
to  define  the  span,  the  greater  is  the  over -est imat i on  of 
stresses  and  displacements.  As  the  number  of  elements  (NSEG) 
is  reduced  from  40  to  10  for  the  span,  there  is  almost  no 
difference  between  the  results.  With  a  6  element  span  the 
errors  are  more  noticeable,  but  still  within  approximately 
10%  deviation  from  the  40  element  span  (these  results  are 
presented  in  Appendix  4  for  Run  Numbers  1  to  4).  This  means 
that  for  the  purposes  of  mine  design  a  6  element  span  is 
acceptable.  For  greater  accuracy  a  10  element  span  will  give 
results  that  are  very  close  to  the  40  element  span  results. 

The  results  are  only  given  for  x  co-ordinates  at  the 
centre  of  each  element,  as  this  is  the  position  for  which 
the  solution  is  derived.  Greater  detail  of  the  stress, 
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strain  and  displacement  distributions  requires  a  larger 
number  of  elements,  to  increase  the  number  of  internal 
points  calculated. 

3.2.3  Simplified  Mine  Layouts 

This  section  contains  the  analysis  of  a  simple  mine 
design  problem.  The  analysis  is  restricted  to  the  variation 
of  three  parameters,  to  facilitate  greater  understanding  of 
the  problem.  In  reality  a  comprehensive  design  approach 
would  include  the  analysis  of  additional  parameters,  and  the 
use  of  other  analysis  methods  for  comparison. 

The  particular  problem  analysed  is  that  of  twin 
headings  (rooms)  driven  parallel  to  each  other  for  the  full 
height  of  a  flat  seam  (Figure  3.3).  The  first  four  program 
runs  (Runs  06,  07,  08  &  09)  investigate  the  effect  of 
variations  in  the  width  of  the  central  pillar  on  the 
stresses  and  displacements  (all  other  parameters  held 
constant).  The  last  two  runs  (Runs  10  <§  11)  investigate  the 
effect  of  variation  of  elastic  parameters  for  the  seam  and 
country  rock,  with  the  geometry  the  same  as  Run  08.  The 
input  data  for  these  runs  can  be  seen  in  Figure  3.4.  The 
horizontal  stress  has  been  assumed  to  be  1/3  of  the  vertical 
stress  by  considering  the  elastic  material  to  be  loaded  by 
the  overburden  with  lateral  constraint  preventing  horizontal 
strain.  The  graphical  output  from  these  program  runs  can  be 
seen  in  Appendix  4  (note  that  metric  units  have  been  used 
after  Run  06).  The  stress  components  for  the  Run  08  results 
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are  presented  in  Figure  3.5,  where  they  may  be  compared 
directly  to  the  physical  layout.  These  diagrams  show  three 
of  the  major  parameters  of  concern;  the  vertical  stress 
within  the  seam,  the  horizontal  stress,  and  the  shear  stress 
within  the  immediate  roof  span. 

Diagram  A  (Figure  3.5)  clearly  shows  the  increase  of 
stress  within  the  central  pillar  and  the  span  abutments.  It 
is  interesting  to  note  that  the  maximum  stress  (pre-failure) 
within  the  pillar  occurrs  at  its  sides,  and  not  at  the 
centre  (this  has  been  demonstrated  by  Wagner,  1974).  The 
average  pillar  stress  is  15.0  MN/m2  greater  than  the  initial 
stress  state  of  25.0  MN/m2,  an  increase  of  the  average 
pillar  vertical  stress  by  60%  of  its  initial  state. 

Diagram  B  (Figure  3.5)  indicates  the  reduction  of 
horizontal  stresses  within  the  immediate  roof  (1.5m  above 
the  top  of  the  seam),  and  the  appearance  of  tensile  stresses 
above  the  pillar.  Tensile  stresses  are  very  undesirable  in 
mine  structures,  as  most  rocks  have  low  tensile  strength. 

Diagram  C  (Figure  3.5)  shows  the  shear  stress  developed 
within  the  immediate  roof  (1.5m  above  the  top  of  the  seam). 
The  maximum  values  occur  above  the  pillar  and  side 
abutments.  These  are  potential  roof  failure  zones 
(especially  with  weak,  thinly  bedded  strata),  where  beam 
type  failures  are  common.  The  combination  of  tensile 
stresses,  together  with  high  shear  stresses  above  the  pillar 
sides,  creates  a  potential  failure  zone.  The  important 
results  from  all  six  runs  are  tabulated  in  Figure  3.6. 
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Analysis  of  the  pillar  and  abutment  stabilities, 
requires  Knowledge  of  the  pillar  strength  characteristics. 
The  vertical  stress  on  the  pillar  can  be  compared  with  its 
estimated  strength  to  calculate  a  safety  factor.  If  an 
acceptable  safety  factor  is  not  achieved,  the  pillar  will  be 
regarded  as  being  unstable.  From  physical  considerations  the 
small  pillar  (width  1.5m)  in  Run  09  would  normally  be 
classified  as  unstable,  because  of  its  low  width  to  height 
ratio  (0.5). 

The  design  of  stable  spans  for  openings  is  often 
undertaken  utilising  past  experience,  and  observation  of 
experimental  openings.  These  analyses  do  not  take  into 
account  the  interactions  with  adjacent  openings,  as  have 
been  investigated  here.  A  simple  analysis  can  be  undertaken 
by  estimating  the  stresses  within  the  immediate  roof  (1.5m 
above  the  top  of  the  seam).  The  horizontal  and  shear 
stresses  at  the  abutments  and  pillar  sides  are  the  major 
factors  for  consideration.  Where  spans  are  large,  the 
stresses  and  displacements  at  the  centre  of  the  span  may  be 
more  important,  as  there  is  no  vertical  support  for  any 
failed  region.  The  results  in  Figure  3.6  indicate  that 
tensile  stresses,  together  with  high  shear  stresses, may 
create  instability  around  the  3m  pillar  (Run  08).  These 
tensile  stresses  are  not  present  in  the  roof  above  the  6m 
pillar  (Run  07)  and  indicates  that  the  6m  pillar  may  give  a 
satisfactory  span  stability. 

Variation  of  the  elastic  moduli  of  the  seam  and  country 
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rock,  were  analysed  to  determine  the  sensitivity  of  the  mine 
structure  (Run  08  -  3m  pillar).  The  effect  on  the  pillar 
stability  was  the  same  for  the  two  Runs  (10  &  11);  doubling 
the  country  rock  stiffness,  or  halving  the  seam  stiffness, 
reduced  the  average  pillar  stress  by  5.4  MN/m2.  The  effect 
on  roof  span  stresses  and  displacements  was  very  small  for 
the  seam  stiffness  variation,  but  the  increase  in  country 
rock  stiffness  reduced  the  displacement  at  the  centre  of  the 
span  by  35%. 

This  simplified  analysis  indicates  that  both  pillar  and 
roof  span  stabilities  are  influenced  by  variations  in  the 
physical  layout.  Variations  in  the  material  constants  have  a 
significant  effect  on  pillar  stability  and  a  marginal  effect 
on  roof  span  stability. 

3.2.4  False  Roof  Tests 

It  has  been  suggested  (Crouch  <&  Fairhurst,  1973)  that 
layered  strata  with  varying  mechanical  properties,  can  be 
modelled  with  unmined  coal  seams.  Because  of  the  nature  of 
the  numerical  model  used,  the  validity  of  this  hypothesis  in 
modelling  false  roof  conditions  was  tested.  A  false  roof 
refers  to  the  existence  of  an  immediate  roof  stratum  with 
different  mechanical  properties  than  the  surrounding  rock. 
The  results  of  Run  16  (identical  to  Run  08,  except  that  a 
false  roof  has  been  modelled  by  another  seam),  can  be  seen 
in  Appendix  4.  Run  16  has  a  false  roof  with  mechanical 
properties  that  are  100  times  stiffer  than  the  country  rock 
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va 1 ues . 

The  particular  simulation  is  of  a  pillar,  delineated  by 
an  opening  on  either  side.  The  creation  of  these  openings 
increases  stresses  within  the  central  pillar.  These  stresses 
should  be  reduced  in  Run  16,  as  a  result  of  the  stiffer  roof 
beam.  The  investigation  of  the  vertical  stresses  at  the 
centre  of  the  pillar  is  adequate,  as  they  are  linearly 
dependent  upon  the  closure  between  roof  and  floor.  Both  sets 
of  results  indicate  a  vertical  stress  at  the  centre  of  the 
pillar  of  39.4  MN/m2 ,  and  hence  the  modelling  of  layered 
strata  using  coal  seams  is  not  a  valid  practise.  Upon 
studying  the  nature  of  the  numerical  model,  these  findings 
can  be  readily  understood.  A  coal  seam  is  modelled  by 
creating  a  hypothetical  crack  that  has  the  top  and  bottom 
attached  to  one  another  by  elastic  spr i nqs .  Upon  excavating 
portions  of  the  seam,  these  spr i nqs  are  removed,  thereby 
removing  the  reactions  between  roof  and  floor.  The  induced 
stress  distribution  is  created  by  the  displacements  which 
take  place  on  these  crack  surfaces.  Therefore  the  presence 
of  a  seam  does  not  change  the  mechanical  properties  of  the 
country  rock,  but  only  the  displacements  which  take  place  on 
the  crack  surfaces.  For  a  seam  which  has  mechanical 
properties  of  infinite  magnitude  one  would  expect  that  no 
displacements  will  take  place  on  the  crack  surfaces,  and 
consequently  there  will  be  no  induced  stresses.  The 
conclusion  is  that  coal  seams  cannot  model  strata  that  are 
stiffer  than  the  country  rock. 
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3.2.5  Multiple  Seam  Problems 

This  section  will  investigate  the  interaction  between 
excavations  and  pillars  in  two  adjacent  seams.  The  results 
demonstrate  the  validity  of  the  numerical  technique  and 
illustrate  simple  interactions. 

Runs  14  &  15  model  two  seams  which  are  relatively  far 
apart,  each  having  two  openings  and  one  central  pillar.  The 
effect  of  offsetting  the  pillar  in  the  lower  seam  is 
investigated  and  the  graphical  results  are  presented  in 
Appendix  4.  The  physical  positions  of  these  two  seams  can  be 
seen  in  Figure  3.7.  The  in-seam  vertical  stresses  are  most 
sensitive  to  the  interaction  mechanism,  and  are  presented  in 
Figure  3.7  for  specific  positions.  It  can  be  seen  that  the 
two  seams  are  sufficiently  far  apart,  in  comparison  to  the 
excavation  size,  so  that  vertical  stress  changes  as  a  result 
of  interactions  are  small. 

The  seams  were  moved  closer  in  Runs  17  &  18  (Results  in 
Appendix  4)  to  increase  the  interaction  effects.  The 
physical  structure,  and  major  vertical  stress  values  for  the 
upper  seam  are  shown  in  Figure  3.8.  The  vertical  stresses 
are  transferred  to  the  right  abutment  in  the  upper  seam,  as 
the  lower  seam  openings  are  moved  to  the  left.  This  results 
from  the  location  of  the  right  excavation  in  the  upper  seam 
above  the  right  abutment  zone  in  the  lower  seam.  This 
increase  in  the  abutment  stress  is  12.5%,  illustrating  the 
adverse  effects  of  offsetting  mine  structures  in  contiguous 
multi-seam  mining. 
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3 . 3  Phvs ica 1  Mode 1 1 i ng 

Problems  associated  with  mine  structures  usually 
involve  some  type  of  rock  failure  mechanism,  that  is 
difficult  to  simulate  with  numerical  modelling  techniques. 

As  in  all  rock  mechanics  problems,  the  most  important 
factors  for  consideration  are  the  rock  mass  parameters,  and 
the  recognition  of  the  failure  mechanisms  involved. 
Sedimentary  rock  structures  are  usually  discontinuous, 
because  of  the  presence  of  bedding  and  jointing.  Physical 
modelling  techniques  have  been  developed  to  simulate  the 
complex  behavioural  mechanisms  which  are  associated  with 
these  strata  types. 

In  physical  modelling  there  are  two  basic  problems;  the 
simulation  of  the  initial  stress  state  by  the  loading 
conditions,  and  the  rock  mass  behaviour.  These  two  problems 
are  interrelated  because  of  the  limitations  of  some  loading 
methods,  especially  those  involving  body  forces  where  the 
use  of  weak  model  materials  is  necessary.  Loading  conditions 
can  be  simulated  by  utilising  friction,  gravity,  momentum  or 
direct  pressure. 

3.3.1  Model  Materials 

The  varieties  of  model  materials  are  numerous 
(Stimpson,  1970)  and  they  seldom  model  the  real  rock  mass 
behaviour  perfectly.  The  material  is  usually  chosen  to  model 
the  parameters  which  are  considered  to  be  of  major 
importance.  Selection  of  scaling  parameters  is  dependent 
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upon  the  predominant  behaviour;  Young's  Modulus  and 
Poisson's  Ratio  for  modelling  elastic  deformations, 
compressive  and  tensile  strengths  for  modelling  plastic 
deformations  (and  failure),  dilatancy  and  friction  angle  for 
modelling  discontinuities. 

3.3.2  Loading  Cond i t ions 

Direct  pressure  is  usually  employed  when  biaxial  or 
triaxial  conditions  are  simulated.  Difficulty  is  associated 
with  producing  boundary  loading  conditions  that  are  capable 
of  differential  displacements  whilst  still  maintaining 
uniform  loading.  This  problem  arises  from  the  physical 
limitations  of  producing  a  model  that  has  the  boundary 
sufficiently  far  from  the  mine  structure  to  avoid 
interact  ions . 

Body  forces  provide  a  means  of  reducing  the  effect  of 
the  boundary  by  loading  each  individual  model  element  in 
proportion  to  its  mass.  Gravitational  loading  is 
proportional  to  density,  which  is  difficult  to  vary,  and 
usually  affects  other  model  parameters.  Centrifugal  loading 
involves  rotation  of  the  model  to  induce  angular 
acceleration  forces  and  provides  the  means  of  varying  the 
loading  conditions  independent  of  the  model  material. 

Impulse  loading  has  been  used  by  vibrating  the  model  in  a 
controlled  manner,  but  this  produces  cyclic  loading 
conditions  and  is  more  suitable  for  earthquake  or  blasting 
vibration  simulation.  Base  friction  loading  is  commonly 
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employed  as  it  allows  more  control  over  the  model  material 
preparation,  and  the  ability  to  stop  and  start  the  model 
test  at  any  point  (Goodman,  1976).  The  body  forces  are 
simulated  by  the  friction  between  the  model  material  and  its 
base.  The  friction  forces  may  be  varied  with  material 
density,  base  friction  factor,  or  by  loading  a  plate  placed 
over  the  mater i a  1 . 

Special  conditions  may  be  simulated  by  combinations  of 
different  loading  methods.  Models  usually  simulate 
two-dimensional  plane  stress  conditions,  because  of  the  ease 
of  preparation  and  observation  of  the  model  response. 


3.3.3  Sea  1 inq  Factors 

In  model  studies,  it  is  not  only  necessary  to  scale  the 
physical  dimensions,  but  also  other  independent  parameters 
involved.  Based  upon  dimensional  analysis,  the  following 
equations  can  be  used  to  match  the  scaling  factors: 


<PL^E)model  “  <PL9/E)actual 


(E/a) 


model 


(E/a) 


actual 


p  =  Density 

L  =  Characteristic  length 
g  =  Gravitational  acceleration 
E  =  Young's  Modulus 
a  =  Stress  conditions 
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Strength  scaling  is  usually  matched  to  the  geometric 
scaling,  or  to  the  loading  conditions.  Dynamic  similitude  is 
not  obtained  in  base  friction  modelling,  but  this  is  not 
usually  considered  to  be  important  in  studying  excavation 
stabi 1 i t i es . 


4.  STABILITY  DESIGN  CONCEPT  ADOPTED 


The  concept  of  utilising  stability  analyses  for  routine  mine 
design  problems  is  presented,  together  with  a  discussion  on 
their  applicability  to  real  mine  situations.  Emphasis  will 
be  placed  on  the  numerical  modelling  technique  utilised 
(DDSEAMS  program);  as  such  techniques  become  more  powerful, 
their  use  as  routine  design  tools  will  increase.  The 
limitations  of  this  approach  will  not  be  overlooked,  but  the 
potential  for  future  development  should  not  be  clouded  with 
excessive  criticism  of  present  capabilities. 


4 . 1  Stress  D istr ibut ions 

Mine  geometries  can  be  seldom  simplified  to  apply 
techniques  such  as  the  tributary  area  theory,  especially 
when  multiple  seam  mining  is  undertaken.  The  numerical 
modelling  technique  adopted  was  used  to  analyse  the 
interactions  induced  between  mine  openings.  The  particular 
program  used  is  a  two-dimensional  formulation,  but  it  could 
be  adapted  for  three-dimensional  analysis.  Elastic  analysis 
is  used,  as  this  allows  a  boundary  element  technique  for 
seam  modelling.  This  reduces  the  input  data  by  one  dimension 
(i.e.  the  2-D  analysis  only  requires  1-D  input  data  on 
opening  geometry).  The  real  behaviour  of  the  rock  mass  is 
usually  unknown,  therefore  it  is  not  necessary  to  undertake 
a  more  complex  analysis. 
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4 . 2  Pillar  Stabi 1 i tv 

The  pillar  stability  is  a  factor  of  the  complete 
stress/strain  (load/displacement)  characteristics  of  the 
pillar.  The  numerical  model  assumes  the  pillars  to  be 
elastic  and  therefore  localised  failures  resulting  in  stress 
re-distributions  are  not  accounted  for.  Hustrulid  (1975) 
presents  pillar  compressive  strength  criteria  based  on 
empirically  derived  size  and  shape  factors  (section  2.1.2). 
The  pillar  strength  is  equivalent  to  the  maximum  average 
stress  that  the  pillar  is  able  to  withstand,  and  therefore 
localised  failures  are  empirically  accounted  for. 

Individual  pillar  stability  can  be  estimated  by 
assuming  that  the  average  stress  on  the  pillar  is  equal  to 
the  average  stress  predicted  by  the  numerical  model.  A 
safety  factor  is  then  defined  as  the  pillar  strength  divided 
by  the  average  pillar  stress.  If  the  safety  factor  is 
greater  than  1.0  then  pillar  failure  is  not  predicted  to 
occur.  Pillar  failure  is  predicted  for  safety  factors  less 
than  1.0,  but  then  the  actual  value  has  no  physical  meaning. 
Stress  re-di str i but  ions  in  the  overall  structure  occur  when 
pillars  fail,  and  these  are  not  accounted  for  in  the 
numer i ca 1  mode  1 . 


4 . 3  Strata  Stress  Cr i ter  ion 

There  are  a  number  of  failure  criteria  that  define  the 
strength  and  failure  mechanism  of  competent  brittle  rock  in 
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triaxial  compression.  The  linear  Mohr /Cou lomb  failure  theory 
will  be  discussed  (Kidybinski  <S  Babcock,  1973),  as  the 
failure  envelope  is  assumed  to  have  a  linear  relationship 
between  shear  and  normal  stresses.  The  stress  state  can  be 
defined  by  drawing  the  Mohr  circle  intersecting  the  normal 
stress  axis  (zero  shear  stress)  at  the  major  and  minor 
principal  stress  values.  An  example  of  a  failure  envelope 
and  a  stress  state  can  be  seen  in  Figure  4.1.  Failure  takes 
place  if  the  Mohr  circle  intersects  or  touches  the  failure 
envelope.  The  definition  of  the  relative  position  of  the 
Mohr  circle,  with  respect  to  the  failure  envelope,  can  be 
used  as  a  basis  for  assessing  the  stability  at  a  point. 

In  attempting  to  define  the  stability  of  a  particular 
stress  state  the  stress  path  to  failure  must  be  defined 
(Kidybinski  <£  Babcock,  1973).  This  path  is  unknown,  and  a 
criterion  that  involves  uniform  expansion  (or  contraction) 
of  the  Mohr  circle  is  employed.  Thus  a  stress  criterion  (SC) 
was  developed  from  Figure  4.1,  given  below: 

SC  =  (Q-R)/Q  =  1  -  (a  —a  3  )  /sine}) .  {o  ^  +  2  .  C/Tan<f>) 

SC  =  STRESS  CRITERION 

Q  =  Distance  from  centre  of  Mohr  circle  to  envelope 
R  =  Radius  of  Mohr  Circle 
a1  =  Major  Principal  stress 
a ^  =  Minor  Principal  stress 
<J>  =  Internal  friction  angle 


C  =  Cohesion 


. 
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This  gives  a  value  of  zero  when  the  point  of  failure  is 
reached,  is  positive  when  conditions  are  stable,  and 
negative  when  failure  has  been  exceeded. 

If  a  safety  factor  is  required  it  can  be  obtained  from: 

Safety  Factor  =  1.0  /  (1.0  -  SC) 

The  use  of  the  stress  criterion  derived  is  adequate  for 
defining  stable  and  unstable  conditions  at  a  point.  A  small 
computer  program  (PSTRESS)  was  developed  to  read  the  results 
from  the  DDSEAMS  program  and  calculate  the  principal 
stresses  and  stress  criterion.  The  program  listing  and 
sample  output  are  presented  in  Appendix  3. 


4 . 4  D i scuss ion  of  Design  Concept  L imi tat  ions 

Some  of  the  analysis  limitations  are  discussed  in  the 
following  chapters,  together  with  comparisons  between  the 
numerical  model,  two  physical  models  and  an  underground 
phenomenological  study.  The  major  limitations  and 
suggestions  for  avoiding  these  are  presented  in  this 
sect i on . 

Critical  analysis  of  the  adopted  approach  requires 
understanding  of  reality,  and  relating  this  to  the  modelling 
simplifications.  If  insufficient  data  is  available  on  the 
rock  mass  parameters  and  post-failure  behaviour,  then 
complex  modelling  using  assumed  values  is  inapplicable.  In 
any  analysis  the  quality  of  the  output  can  never  be  better 
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than  the  quality  of  the  input. 

Previous  work  (Ortlepp  &  Cook  ,1964;  Ortlepp  &  Nicoll, 
1964  and  Salamon  &  Oravez,  1970)  has  indicated  that 
stratified  rock  masses  can  be  simulated  as  linear  elastic 
bodies  with  reasonable  accuracy.  Starfield  &  McClain  (1973) 
have  investigated  non-linear  behaviour  (creep)  in  salt 
pillars,  loaded  by  an  elastic  rock  mass,  and  achieved 
acceptable  results.  These  results  indicate  that  if, 
non-linear  behaviour  in  the  proximity  of  the  mine  structure 
is  modelled,  the  remainder  of  the  rock  mass  can  be  assumed 
to  be  linearly  elastic.  The  remaining  problem  is  to  define 
the  non-linear  and  post-failure  behaviour  of  the  rock  mass 
in  the  vicinity  of  the  mine  structure.  The  modelling  of 
these  areas  is  more  suited  to  finite  element  techniques 
(Kidybinski  &  Babcock,  1973),  as  different  parameters  can  be 
allocated  to  regions  according  to  the  pre-defined  degrees  of 
failure.  The  boundary  element  program  formulation  may  be 
adapted  to  model  non-linear  behaviour  within  the  seam,  but 
the  country  rock  must  have  a  linear  elastic  behaviour. 

The  advantage  of  using  the  DDSEAMS  program  is  that  the 
sensitivity  of  the  mine  structure  to  parametric  and 
geometric  variations  can  be  investigated  relatively  easily. 
The  choice  of  a  preferred  geometry  can  be  made,  on  the  basis 
of  defined  factors  (e.g.  minimising  unstable  regions),  by 
analysing  various  mine  geometries.  If  the  analysis  indicates 
that  failure  will  occur,  then  the  results  will  not  be  a  true 
model  of  reality.  The  occurrence  of  failures  involve  stress 
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redistributions  that  are  not  accounted  for  in  the  numerical 
model.  A  more  accurate  prediction  can  be  obtained  by  finite 
element  analysis  or  physical  modelling  with  correct 
s i mi  1 i tude . 

Back  analysis  techniques  should  be  used  wherever 
possible  to  ensure  that  the  model  adequately  simulates  the 
real  rock  mass.  Where  time  dependent  behaviour  is 
encountered  a  further  dimension  is  added  to  the  problem, 
greatly  increasing  the  complexity. 

The  strata  stress  criterion  only  indicates  potential 
failure  zones,  and  not  the  stability  of  the  overall 
structure.  Therefore  it  does  not  define  the  final  failed 
regions,  as  stress  re*di s t r i but i ons  as  a  result  of  failure 
are  not  accounted  for . 

Pillar  safety  factors  are  adequate  for  stable 
conditions  (all  pillar  safety  fa'ctors  greater  than  1.0).  If 
any  pillars  have  safety  factors  less  than  1.0,  the  values 
for  the  remainder  of  the  pillars  do  not  account  for  the 
stress  re-distributions  that  will  occur. 


. 
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5.  MULTI -SEAM  STUDIES 


A  comparative  study  was  undertaken  to  model  more  complex 
multi-seam  interactions  with  the  numerical  technique  and 
compare  the  results  with  two  physical  models  and  an 
underground  phenomenological  study.  No  suitable  mine  data 
was  available,  consequently  the  simulation  presented  was 
designed  primarily  for  physical  and  numerical  model 
comparisons.  The  layout  presented  is  based  upon  an  operating 
mine  at  Grande  Cache,  Alberta.  The  upper  seam  position  has 
been  moved  closer  to  the  lower  seam,  to  match  the  dimensions 
of  the  physical  model  apparatus,  and  to  increase  the  seam 
interactions.  Therefore  the  model  geometry  is  not  a  perfect 
representat i on  of  the  actual  mine  situation. 


5 . 1  Numer ica 1  Model 

The  geometric  layout  for  the  model  is  presented  in 
Figure  5.1,  where  large  openings  have  been  excavated  in  the 
lower  seam  to  simulate  de-pi  1  Taring  operations.  Openings  in 
the  upper  seam  have  been  introduced  to  simulate  development 
work  concurrent  with  undermining  operations.  The  results  for 
the  two  seams  are  presented  graphically  in  Figure  5.2  and 
Figure  5.3,  together  with  the  input  parameters. 

A  grid  system  of  data  points  was  set  up,  and  the 
analysis  results  from  these  points  were  used  to  develop  a 
series  of  contour  plots  of  major  results.  These  contour 
plots  only  depict  the  results  for  the  country  rock.  The 
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vertical  and  horizontal  stress  contours  are  presented  in 
Figure  5.4  and  Figure  5.5  respectively.  The  analysis  results 
were  used  as  input  data  to  the  PSTRESS  program  and  principal 
stress,  maximum  shear  stress  and  stress  criterion  values 
were  produced.  Contour  plots  of  major  principal  stress, 
minor  principal  stress,  maximum  shear  stress  and  stress 
criterion  values  are  presented  in  Figure  5.6,  Figure  5.7, 
Figure  5.8  and  Figure  5.9  respectively.  A  contouring  package 
available  at  the  University  of  Alberta  (Surface  II)  was  used 
to  produce  these  plots.  The  smoothness  of  these  plots  is  a 
function  of  the  coarse  grid  pattern  adopted.  The  potentially 
unstable  areas  predicted  by  negative  stress  criterion  values 
are  shaded  in  Figure  5.9. 

The  vertical  stress  contours  (Figure  5.4)  indicate  the 
induced  stresses  as  a  result  of  excavating  the  mine 
openings.  The  initial  vertical  stresses  at  the  upper  and 
lower  seam  depths  are  11.5  and  12.5  MN/m2  respectively.  The 
high  stress  gradients  at  the  pillar  abutments  are  evident 
(the  smoothness  of  the  stress  contours  in  these  regions  is 
reduced  because  of  the  contour  package  used).  The  vertical 
stress  relief  above  the  openings  is  more  pronounced  above 
opening  No. 4  in  the  lower  seam,  as  a  result  of  opening  No. 1 
located  above  it  in  the  upper  seam.  The  average  vertical 
stress  on  the  pillars  (tabulated  on  page  48)  indicates  that 
the  small  pillar  between  openings  3  and  4  has  an  average 
vertical  stress  of  33.0  MN/m2.  If  this  average  stress  had 
been  estimated  by  the  tributary  area  theory,  then  an  average 
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stress  of  60.4  MN/m2  would  be  expected.  This  indicates  the 
improved  stress  prediction  capabilities  of  the  numerical 
model  which  accounts  for  the  proximity  of  the  abutments. 

The  horizontal  stress  contours  (Figure  5.5)  also 
indicate  the  stress  relief  above  the  openings,  but  to  a 
lesser  degree  than  the  vertical  stress  relief.  High 
horizontal  stress  gradients  below  opening  No.1,  combined 
with  low  vertical  stresses,  indicate  potential  instability 
that  is  confirmed  in  the  stress  criterion  contour  plot. 

The  major  principal  stress  contours  (Figure  5.6)  and 
minor  principal  stress  contours  (Figure  5.7)  indicate  that 
the  horizontal  stress  is  the  major  principal  stress  in  zones 
of  vertical  stress  relief.  In  zones  of  high  vertical  stress 
the  vertical  stress  is  the  major  principal  stress. 

The  maximum  shear  stress  contours  (Figure  5.8)  indicate 
the  high  shear  stress  gradients  at  the  corners  of  openings. 
These  are  also  predicted  by  the  stress  criterion  to  be  areas 
of  potential  failure.  The  stress  criterion  is  based  on 
intact  rock  shear  failure  and  therefore  the  maximum  shear 
stress  is  the  major  variable  which  indicates  possible 
failure  zones.  However,  it  does  not  allow  for  the  influence 
of  the  principal  stresses  and  therefore  does  not  give  a 
complete  indication  of  the  potential  for  failure. 

The  stress  criterion  contours  (Figure  5.9)  gives  a  much 
improved  prediction  of  over -stressed  regions.  The 
over -stressed  region  below  opening  No.1  is  evident,  and  this 
region  has  not  been  predicted  to  be  over-stressed  by  the 
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maximum  shear  stress  alone.  The  small  over-stressed  zones  at 
the  extremities  of  the  upper  seam  are  a  result  of  boundary 
condition  limitations  of  the  numerical  model. 

The  material  strength  parameters  have  been  simplified 


to  give  zero 

cohesion  and 

i nterna 1 

friction  angles  given 

below: 

Mater ia 1 

Locat ion 

Internal  Friction  Anale 

MUDSTONE 

above  top  seam 

30.0° 

SILTSTONE 

between  seams 

32.0° 

SANDSTONE 

below  lower 

seam 

35.0° 

These  assumptions  were  made  on  the  basis  of  published  data 
( Hoek  <S  Brown,  1980)  for  these  rock  types,  as  no  test  data 
was  available.  The  internal  friction  angle  was  taken  from 
the  tangent  to  the  failure  envelope  at  a  confining  stress 
equal  to  50%  of  the  uniaxial  compressive  strength. 


5 . 2  Base  Fr ict ion  Model 

One  base  friction  model  run  was  made,  primarily  to 
investigate  the  caving  char acter i st i cs  of  bedded  strata.  The 
base  friction  apparatus  consists  of  an  electrically  driven 
sand  paper  belt  (Plate  1).  A  flour  and  methanol  model 
material  was  used,  that  was  easy  to  prepare  and  cut.  The 
layout  prior  to  excavation  (Figure  5.10)  is  similar  to  the 
numerical  model,  except  that  an  additional  opening  was  made 
in  the  upper  seam.  The  model  material  is  stronger  than  the 
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stresses  imposed  by  the  base  friction  loading,  and  therefore 
failure  along  discontinuities  is  predominant.  The 
horizontally  bedded  strata  were  modelled  by  cutting  beds  at 
5.1cm  intervals,  with  two  3.8cm  beds  above  the  lower  seam, 
and  a  seam  thickness  of  5.1cm.  The  real  height  of  the 
opening  is  3.5m,  representing  a  geometric  scale  factor  of 
1/69.  The  development  openings  in  the  upper  seam  were 
excavated  first,  followed  by  the  two  openings  on  the  right 
hand  side  of  the  lower  seam.  The  main  openings  (No. 3  and 
No . 4 )  in  the  lower  seam  were  then  excavated  with  no 
significant  displacements  till  the  openings  reached  their 
full  spans  (Plate  4).  The  initial  roof  failure  occurred 
above  these  spans;  with  the  two  immediate  beds  failing  and 
significant  bed  separation  of  the  two  above.  The  small 
central  pillar  in  the  lower  seam  was  then  excavated,  to 
simulate  failure,  and  the  resultant  cave  is  shown  in 
Plate  5.  Significant  bed  separation  up  to  six  beds  above  the 
lower  seam  was  observed,  with  the  subsidence  profile  in  the 
upper  strata  well  pronounced.  A  cave  angle  of  approximately 
70°  was  observed  at  both  abutments. 

The  base  friction  model  is  very  qualitative  and  cannot 
be  used  to  accurately  predict  the  real  situation.  However, 
it  does  serve  as  a  useful  model  to  indicate  the  behaviour  of 
discontinuous  strata  around  excavations. 


. 
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5 . 3  3-D  Phvs ica 1  Model 

The  3-D  physical  model  was  developed  to  simulate  the 
fracture  char acter i st i cs  of  pillars  (SzwilsKi  &  Whittaker, 
1975)  and  to  observe  multi-seam  interactions  qualitatively. 
The  coal  pillars  are  modelled  by  p 1  aster -of -par i s  cement  and 
the  rock  strata  by  a  series  of  rubber  strips.  The  rubber 
strips  offer  a  lateral  strain  effect  and  a  high 
rubber/pi  1 lar  ( strata/pi  1 1 ar )  contact  friction  factor.  The 
model  is  loaded  uniaxial ly  in  a  200  ton  compression  frame 
(Plate  3 )  by  a  hydraulic  ram  connected  to  a  braced  steel 
plate  0.91m  by  0.61m. 

Two  models  were  constructed,  Run  1  had  5.1cm  thick 
seams,  and  run  2  had  3.8cm  thick  seams. The  horizontal  scale 
factor  was  1/120  for  both  runs  and  the  vertical  scale 
factors  were  1/60  for  Run  1  and  1/90  for  Run  2.  The  plaster 
mix  used  had  a  uniaxial  compressive  strength  of  5.2  MN/m2. 
Attempts  were  made  to  produce  weaker  pillars  of  plaster/sand 
mixes,  but  these  were  abandoned  after  encountering  problems 
with  pouring  and  handling  large  pillars.  The  pillar  layouts 
in  the  two  seams  are  shown  in  Figure  5.11  (both  runs 
identical)  and  the  assembly  for  the  Run  2  model  is  shown  in 
Plates  9,  10  &  1  1  . 

The  model  was  progress i ve 1 y  loaded  and  photographs  were 
taken  at  intervals,  until  a  yield  state  was  reached.  The 
yield  state  was  indicated  by  the  peak  load  condition, 
associated  with  increased  displacement  of  the  model.  The 
yield  state  of  Run  1  was  at  20  MN/m2  ram  pressure 
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(0.96  MN/m2  average  model  stress),  shown  in  Plate  6.  The 
deformation  of  opening  No. 2  is  shown  in  Plate  7;  floor 
heave,  roof  sag  and  pillar  slabbing  failures  were  observed. 

After  failure  the  model  was  dismantled  and  the  pillar 
failures  were  observed  in  detail.  A  plan  view  of  openings 
No. 1  and  No. 2  in  the  upper  seam  is  shown  in  Plate  8.  Note 
that  there  has  been  a  much  greater  slabbing  failure  of 
pillar  sides  in  opening  No . 2 . 

Run  2  was  loaded  in  a  similar  manner  and  pictures  of 
model  average  stresses  of  0.78  MN/m2  and  0.92  MN/m2  (yield 
state)  are  shown  in  Plates  12  <5  13  respectively.  One  of  the 
pillars  between  openings  No. 5  and  No. 6  in  the  lower  seam, 
shown  in  Plate  14,  illustrates  the  pillar  fractures  running 
parallel  to  the  pillar  sides.  This  pillar  failure  mode  is 
further  emphasized  in  the  upper  seam  between  openings  No. 1 
and  No . 2  (Plate  15). 


5 . 4  Underground  Phenomeno logical  Study 

The  subject  of  this  study  is  a  coal  mine  located  at 
Grande  Cache,  Alberta.  The  coal  measures  are  Lower 
Cretacious,  occur  in  the  Rocky  Mountains,  and  are  situated 
in  a  flat-bottomed  syncline.  The  inter-seam  strata  consists 
of  interbedded  sandstones,  siltstones  and  shales.  Operations 
are  carried  out  in  the  No. 4  seam  (lower)  and  the  No. 11  seam 
(upper).  The  intermediate  seam  (No. 10),  simulated  in  the 
physical  and  numerical  models,  has  not  been  represented 
correctly  in  the  vertical  dimension. 
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Mining  is  carried  out  by  room  and  pillar,  coal  being 
extracted  by  continuous  miners.  Panels  are  developed  on  a 
three  entry  system,  rooms  being  developed  on  advance  and 
pillars  extracted  on  retreat.  Rooms  are  6m  wide  and  3m  high, 
with  the  roof  supported  by  bolting  on  a  1.2m  square  pattern. 
Rib-sides  are  supported  by  rib  posts  and  lagging.  Pillar 
width  depends  upon  conditions,  but  is  usually  24.4m  (30.4m 
roadway  centres). 

The  initial  development  in  the  No. 4  seam  (6.1m  thick) 
is  confined  to  the  upper  3m.  However,  the  lower  portion  of 
the  seam  is  removed  during  de-pi  1  Taring  operations,  and  the 
roof  is  allowed  to  cave.  Regular  roof  caving  is  promoted  by 
avoiding  remnant  pillars,  and  maintaining  a  uniform  break 
line  for  the  cave.  Irregular  caving  is  associated  with 
excessive  pillar  stresses,  floor  heave  and  instability  of 
roadway  intersections  close  to  the  break  line. 

A  newly  exposed  rib-side,  pictured  in  Plate  16,  shows 
the  rib  posts,  lagging  and  bolting  with  wire  mesh.  The 
competence  of  the  pillar  is  noticeable,  together  with  the 
minor  seam  irregularities  indicated  by  the  roof  "roll".  The 
initial  indications  of  high  pillar  stresses  are  seen  by  the 

pillar  expansion  along  a  soft  band  close  to  the  roof 

(Plate  17).  The  next  stage  in  pillar  failure  is  indicated  by 
fractures  opening  parallel  to  the  pillar  sides  (Plate  18), 
causing  pillar  expansion  (Plate  19).  Pillar  expansion  and 
floor  heave  lead  to  failure  of  rib  supports  (Plate  20).  The 

floor  heave  creates  closure  of  the  openings  and,  in  order  to 
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maintain  the  working  height,  the  floor  has  to  be  excavated 
(brushing).  Floor  heave  close  to  the  cave  line  is  indicated 
in  Plate  21.  Regular  caving  to  the  break  line  reduces 
pressure  on  supports,  as  is  shown  in  Plate  21.  The  thin 
layered  roof,  shown  in  Plate  22,  indicates  a  cave  of  the 
immediate  roof  at  an  intersection.  It  is  interesting  to 
observe  the  failure  initiated  between  roof  bolt  anchor 
plates,  with  few  anchor  failures.  The  development  opening  in 
the  upper  portion  of  the  seam,  together  with  the  mined  lower 
portion,  is  shown. 

It  is  clear  from  this  study  that  the  real  mine 
situation  exhibits  time-dependent  pillar  and  floor  failures 
in  highly  stressed  regions. 


5 . 5  Comparison  of  results 

The  vertical  stress  contour  plot  (Figure  5.4) 
illustrates  the  re-distribution  of  vertical  stresses  around 
the  mine  openings.  Stress  concentrations  are  induced  within 
the  pillars,  especially  those  between  openings  No. 3,  No. 4, 
No. 5  and  No. 6  in  the  lower  seam.  The  vertical  stress  relief 
above  the  openings  is  observed,  especially  above  opening 
No. 4,  as  a  result  of  the  presence  of  another  opening  above 
it.  The  vertical  stress  on  the  coal  pillars  is  taken  to  be 
equivalent  to  the  average  vertical  stress  (predicted  by  the 
numerical  model),  uniformly  distributed  over  the  pillar.  The 
pillar  strength  can  be  estimated  (Jeremic,  1979)  from  the 
following  equations: 
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a  =  c.W/H  (W/H  <  5) 

P 

a  =  2.4  +  c.W/H  (W/H  >  5) 

ir 

a  =  Pillar  compressive  strength  (MN/m2 ) 

W  =  Width  of  pillar 

H  =  Height  of  pillar 

c  =  Factor  depending  upon  coal 

(varies  from  1.2  to  1.6  MN/m2) 

A  value  of  1.5  MN/m2  was  assumed  for  C  in  all  calculations. 
Pillar  widths  are  assumed  to  be  24m  along  the  long  axis  of 
the  openings,  from  which  the  effective  width  can  be 
calculated: 

W  =  /24 . W “ 
e  t 

W^  =  Effective  pillar  width 

W  =  Width  of  pillar  transverse  to  long  axis 

The  following  table  gives  the  average  pillar  vertical 
stress,  effective  width,  width  to  height  rat i o , s trength ,  and 
safety  factor  (the  pillars  being  defined  with  respect  to  the 
openings  e.g.  pillar  1/2  is  the  pillar  between  openings  No. 1 
and  No . 2 .  )  : 
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Pillar 

Vert ica 1 

Stress 

(MN/m2 ) 

Effect i ve 

Width 

(m) 

W/H 

Strenath 

(MN/m2 ) 

Safetv 

Factor 

1/2 

12.5 

20.8 

6.9 

12.8 

( 1.02)* 

3/4 

33.0 

12.0 

4.0 

6.0 

(0.18)* 

4/5 

19.4 

24.0 

8.0 

14.4 

(0.74)* 

5/6 

21.2 

14.7 

4.9 

7.4 

(0.35)* 

*  The  results  indicate  that  progressive  failures, 
resulting  in  stress  re-distributions,  will  occur.  Therefore 
the  safety  factors  calculated  do  not  represent  final  values, 
but  only  predict  instability.  It  should  be  noted  that  the 
instability  in  the  lower  seam  is  associated  with  a  caving 
phenomenon,  that  occurs  on  retreat,  and  therefore  safety 
factors  below  1.0  would  be  expected. 

It  was  observed  that  opening  No. 2  of  the  3-D  physical 
model  Run  1  (Plates  657)  had  characteristic  roof  and  floor 
displacements.  For  a  single  opening  in  elastic  rock,  it 
would  be  expected  that  the  roof  and  floor  displacements, 
causing  closure,  would  be  equal  and  symmetrical.  Upon 
inspection  of  Plate  7  it  was  observed  that  the  floor  heave 
was  much  greater  than  the  roof  sag.  The  numerical  model  has 
much  stiffer  elastic  constants  than  the  rubber ,  but  a 
comparison  between  the  two  was  made  by  applying  horizontal 
and  vertical  scale  factors.  The  numerical  model 
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displacements  were  scaled  to  match  the  outline  of  the 
physical  model  (Figure  5.12).  This  qualitative  comparison 
indicates  that  the  numerical  model  also  predicts  that  the 
floor  heave  is  greater  than  the  roof  sag,  but  to  a  lesser 
degree  than  the  physical  model.  Close  inspection  of  Plate  6 
indicates  bed  separation  occurring  in  the  floor  of  the  model 
opening,  accounting  for  the  increase  in  floor  displacements. 

The  contour  diagram  of  stress  criterion  (Figure  5.9)  is 
most  informative,  as  it  combines  all  the  data  on  the  major 
principal  stress  (Figure  5.6),  minor  principal  stress 
(Figure  5.7)  and  maximum  shear  stress  (Figure  5.8),  together 
with  the  strength  properties  of  the  material.  The 
instability  of  opening  No. 2  is  clearly  indicated  by  the 
shaded  zones  in  Figure  5.9,  and  the  physical  model 
(Plates  6  &  8 )  demonstrates  the  pillar  slabbing  along  this 
opening . 

The  mechanism  of  pillar  failures  by  fracturing  parallel 
to  pillar  sides  is  shown  in  Plates  14  <£  15,  with  the 
constrained  pillar  core  providing  most  of  the  support.  The 
underground  photographs  show  a  similar  failure  mechanism 
(Plates  18  &  20) . 

Opening  stability  is  primarily  concerned  with  roof  span 
failures,  with  roof  support  by  rock  bolting  being  commonly 
used  in  room  and  pillar  mining.  The  existence  of  layered 
strata  in  the  immediate  roof  of  the  underground  mine 
(Plate  22),  indicates  that  the  analysis  of  stability  should 
be  undertaken  assuming  that  bolting  forms  a  composite  roof 
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beam.  It  is  interesting  to  note  that  the  majority  of  the 
roof  cave  depicted  (Plate  22)  has  been  caused  by  rock 
failure  between  bolt  anchor  plates.  If  roof  support  produces 
very  competent  roof  beams  it  may  have  a  detrimental  effect 
on  the  mine  conditions  during  full  extraction,  as  regular 
caving  may  be  prevented. 


6.  MINE  DESIGN  SIMULATION 


A  mine  design  simulation  was  undertaken  using  data  from 
proposed  underground  mining  of  the  Silkstone  seam  at  Coal 
Valley,  Alberta.  The  study  is  intended  to  illustrate  the 
application  of  the  mine  design  concept  adopted  and  is  too 
limited  in  information  and  scope  to  be  regarded  as  a 
comprehensive  stability  analysis. 

Information  on  structural  geology  was  obtained  from  the 
mine  personnel  and  rock  mass  parameters  were  estimated  from 
a  geotechnical  testing  program  on  drill  core  from  the 
proposed  site. 


6 . 1  Geo  1 oqy 

Only  a  brief  summary  will  be  presented  here  as  the 
limited  nature  of  this  study  does  not  warrant  a  detailed 
descr i pt i on . 

6.1.1  Structure 

In  the  proposed  area  the  predominant  feature  is  a  broad 
syncline  with  dips  varying  from  4°  to  18°,  with  an  average 
of  10°,  towards  the  southwest.  The  area  is  bounded  by  major 
fault  systems  with  faulting  dividing  the  area  into  five 
offset  blocks.  The  Silkstone  seam  outcrops  along  the 
northeast  limb  of  this  syncline,  which  plunges  towards  the 
southeast.  Previous  drilling  has  not  indicated  the 
predominance  of  minor  faulting  and  seam  "rolls",  but  this 
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may  be  due  to  the  low  density  of  drilling.  The  area  is 
bounded  on  the  down-dip  side  by  an  undefined  fault  (approx 
1500m  from  the  outcrop),  where  the  seam  is  at  a  depth  of 
approximately  150m. 

6.1.2  Strat iqraohv 

The  Si IKstone  coal  measures  consist  of  two  seams,  the 
upper  denoted  as  the  Wee  and  the  lower  as  the  Bourne, 
separated  by  approximately  9m  of  silty  mudstones.  The  Bourne 
seam  averages  only  0.9m  in  thickness  and  is  not  considered 
to  be  economically  recoverable.  The  average  normal  thickness 
of  the  Wee  seam  is  3m,  but  varies  from  1.5m  to  5.8m.  In 
places  the  seam  is  divided  by  a  carbonaceous  shale  and 
mudstone  parting  which  sometimes  exceeds  0.6m.  The  strata 
above  the  Wee  seam  consists  mainly  of  massive  sandstone  and 
siltstone  beds.  In  the  eastern  portion,  close  to  the  outcrop 
(where  three  of  the  four  core  holes  were  located),  the  roof 
strata  is  predominantly  interbedded  siltstones  and 
mudstones . 

The  coal  is  ranked  as  High  Volatile  Bituminous  C  by  the 
ASTM  classification. 


6 . 2  Geotechnical  Test ing  Program 

Geotechnical  testing  of  core  from  four  drill  holes  was 
undertaken  to  define  the  rock  mass  parameters  and  other 
properties  which  relate  to  the  mine  environment.  The  core 
was  obtained  for  shallow  seam  conditions,  as  it  was  proposed 
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to  initially  undertake  underground  exploration.  Conditions 
at  greater  depth  may  be  different,  but  no  core  was  available 
for  testing. 

The  diamond  drilling  for  2.6  inch  diameter  core  was 
carried  out  using  water  flush  and  a  triple  barrel.  Runs  were 
of  approximately  7  ft.  and  the  core  was  immediately  sealed 
in  the  plastic  inner  tubing.  The  sealed  core  tubes  were 
transported  to  the  laboratory  where  geotechnical  logging  and 
sampling  was  carried  out.  Upon  exposure  to  the  atmosphere 
the  samples  were  sealed  and  placed  in  a  humid  environment 
until  they  were  opened  for  preparation  and  testing.  The 
geotechnical  core  logs  for  the  four  holes  can  be  seen  in 
Appendix  5.  A  detailed  description  of  the  testing  procedure 
is  presented  in  Appendix  6,  and  the  detailed  testing  results 
are  presented  in  Appendix  7.  A  summary  of  the  test  results 
fol lows : 

6.2.1  Uniaxial  Compression 

Samples  of  the  full  core  diameter  were  tested  in 
uniaxial  compression  (with  strain  gauges  attached),  however, 
no  suitable  coal  samples  were  available  for  this  test.  The 
results  obtained  were  used  to  estimate  the  Young's  Modulus, 
Poisson's  Ratio  and  uniaxial  compressive  strength.  Six 
samples  were  tested  and  the  average  results  are  tabulated 
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Youna' s 

Poisson' s 

Compressive 

Modulus 

(MN/m2 ) 

Rat  io 

Strenath 

(MN/m2 ) 

SANDSTONE 

1.36  x  104 

0.267 

73.9 

(3  samples) 

SILTSTONE 

1.00  x  104 

0.288 

66.2 

(2  samples) 

MUDSTONE 

8.62  x  103 

0.172 

70.2 

(1  sample) 

COAL  * 

1.92  x  103 

- 

15.6 

(1  sample) 

*  under cored  sample  (no  strain  gauges  attached) 

6.2.2  Brazi 1 ian  Disc 

Disc  samples  (thickness  approximately  half  of  the 
diameter)  of  both  1.6  inch  and  2.6  inch  diameter  were  used 
in  the  indirect  tensile  strength  Brazilian  tests.  The 
average  results  from  35  tests  are  tabulated  below: 


SANDSTONE  (14  samples)  =  2.29  MN/m2 
SILTSTONE  (6  samples)  =  3.47  MN/m2 
MUDSTONE  (11  samples)  =  3.97  MN/m2 
COAL  (4  samples)  =  0.91  MN/m2 
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6.2.3  Tr iaxia 1  Compression 

In  order  to  investigate  the  compressive  strength  under 
confined  conditions,  and  to  develop  a  Mohr/Coulomb  failure 
envelope,  21  samples  were  tested  in  triaxial  compression. 
Core  samples  were  prepared  by  undercoring  to  1.6  inch 
diameter  to  fit  the  triaxial  cell,  and  the  Jjn  situ  moisture 
content  was  preserved  as  much. as  possible.  There  were 
insuffient  samples  to  produce  a  failure  envelope  for 
siltstone,  but  the  Mohr's  circles  and  failure  envelopes  for 
sandstone,  mudstone  and  coal  were  drawn  (Figures  6.1,  6.2 
and  6.3).  It  should  be  noted  that  some  of  the  height  to 
width  ratios  tested  were  below  2,  because  of  sample 
preparation  problems,  but  the  results  from  these  samples 
appear  to  be  reasonable. 

From  the  results,  simplified  linear  Mohr/Coulomb 
failure  envelopes  have  been  developed,  which  have  the 
following  relationships  for  confining  stresses  less  than 
15  MN/m2: 


O 


SANDSTONE  : 


x  =  6.9  +  cfTan45 


t  =  6.9  +  aTan38 


O 


MUDSTONE : 


o 


COAL: 


t  =  3.0  +  crTan45 


t  =  Shear  stress  (MN/m2 ) 


a  =  Normal  stress  (MN/m2 ) 


' 
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The  low  values  for  coal  cohesion  can  be  explained  by  the 
presence  of  cleat  planes  which  results  in  a  discontinuous 
material.  The  influence  of  bedding  angle  was  not 
investigated  as  all  samples  were  undercored  along  the 
original  core  axis,  hence  the  principal  stress  direction  was 
approximately  80°  to  the  bedding  dip. 

6.2.4  D i rect  Shear 

Shear  strength  along  discontinuities  was  investigated 
using  a  direct  shear  apparatus.  Only  a  small  number  of  tests 
were  carried  out,  because  of  a  shortage  of  suitable  samples. 
Two  mudstone  joints,  one  mudstone  and  one  coal  bedding  plane 
were  sheared,  all  shearing  taking  place  along  the  dip  of  the 
discontinuity.  There  were  no  significant  peak  values  and  the 
residual  values  are  given  below: 

MUDSTONE  JOINT  (50°  dip):  Friction  angle  =  19.9° 

MUDSTONE  JOINT  (30°  dip):  Friction  angle  =  27.1° 

MUDSTONE  BEDDING  (15°  dip):  Friction  angle  =  23.3° 

COAL  BEDDING  (15°  dip):  Friction  angle  =  23.8° 

The  low  friction  angles  in  comparison  to  the  internal 
friction  angles  are  a  result  of  testing  discontinuity 
weakness  planes. 

6.2.5  Slake  Durabi 1 i tv 


Slake  durability  tests  on  the  immediate  floor  strata 


. 
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were  undertaken  to  assess  the  effect  of  water  on  the 
competence  of  these  rock  units.  From  the  standard  testing 
procedure,  results  of  the  Slake  Durability  Index  (second 
cycle)  were  16.8%,  78.5%  and  95.3%  for  three  footwall 
mudstone  samples.  The  higher  values  appeared  to  be 
associated  with  the  sand  content  within  the  sample  and  hence 
the  consequent  reduction  in  clay  content. 

6.2.6  Swel 1 inq  Index 

Two  Swelling  Strain  Index  (unconfined)  tests  were 
carried  out,  one  floor  mudstone  and  one  roof  siltstone.  The 
results  were  2.9%  and  2.5%  for  the  siltstone  and  mudstone 
respectively.  Partial  disintegration  of  the  siltstone  and 
complete  disintegration  of  the  mudstone  took  place,  this 
indicates  that  the  maximum  swelling  would  probably  have  been 
greater  than  that  recorded. 


6 . 3  Rock  Mass  Parameters 

The  rock  mass  parameters  that  are  required  for  analysis 
purposes  are;  the  relationship  between  stress  and  strain 
(Young's  Modulus  <£  Poisson's  Ratio),  and  the  strength 
characteristics.  Estimates  of  these  parameters  have  been 
made  from  laboratory  testing  of  intact  samples,  which  tend 
to  be  stiffer  and  stronger  than  the  rock  mass  itself. 

Scaling  factors  need  to  be  applied  from  the  laboratory 
results  to  obtain  j_n  s i  tu  estimates.  This  process  is  not 
well  understood  as  the  real  rock  mass  is  usually 
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non -homogeneous  and  discontinuous,  thus  forming  a  complex 
structure.  There  is  insufficient  information  to  undertake  a 
detailed  analysis,  therefore  simple  scaling  factors  will  be 
app lied. 

Based  on  previous  work  (Raney  e_t  aj_,  1976)  the  Young's 
Modulus  for  rock  and  coal  was  reduced  by  factors  of  3  and  5 
respectively.  Poisson's  ratio  estimates  from  laboratory 
testing  are  not  regarded  as  being  representative  of  the  rock 
mass,  hence  values  of  0.19  and  0.3  for  rock  and  coal  will  be 
used  (Pariseau  &  Sorensen,  1979  and  Kulhawy,  1975).  The 
Shear  Modulus  will  be  estimated  from  the  Young's  Modulus  and 
Poisson's  Ratio  (assuming  an  isotropic  elastic  material), 
using  the  following  formula  (Plane  Strain  analysis): 

G  =  E/2  (1+v) 

G  =  Shear  Modulus 
E  =  Young ' s  Modulus 
V  =  Poisson's  Ratio 


The  rock 

previous 

mass  parameters 

assumptions  are 

from  the  laboratory 

given  below: 

tests  and  the 

Material 

E  (MN/m2) 

G  (MN/m2) 

ROE 

ROCK 

3.9  x  103 

1.6  x  103 

0.19 

COAL 

3.9  x  102 

1.5  x  102 

0.30 
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E  =  Young's  Modulus 
G  =  Shear  Modulus 
ROE  =  Poisson's  Ratio 

From  triaxial  testing  the  linear  Mohr/Coulomb  failure 
envelopes  have  been  estimated.  These  estimates  do  not  take 
into  account  the  effect  of  discontinuities  in  the  rock  mass, 
and  hence  they  have  been  adjusted  to  give  the  following 
re  1  at i onsh i ps : 


Material 

C  ( MN/m2 ) * 

PHI  (dearees)** 

SANDSTONE 

2.3 

40.0 

MUDSTONE 

2.3 

33.0 

COAL 

0.55 

40.0 

C  =  Rock  Mass  Cohesion 

PHI  =  Internal  Friction  Angle 

*  The  cohesion  values  were  obtained  by  dividing  the  rock 
and  coal  cohesion  values  (from  testing)  by  3  and  5 
respect i ve 1 y . 

**  The  internal  friction  angles  were  obtained  by 
subtracting  5°  from  the  laboratory  results  to  allow  for  the 
low  discontinuity  friction  angles  obtained  from  the  direct 


, 
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shear  tests. 

These  adjustments  also  bring  the  estimates  closer  to  values 
used  in  previous  published  work  (Kidybinski  &  Babcock,  1973 
and  Fadeev  &  Abdyldayev,  1979). 

The  effects  of  anisotropy  and  time  dependent  behaviour 
have  not  been  considered  as  they  are  difficult  to 
investigate  and  the  required  complexity  is  beyond  the  scope 
of  this  study. 


6 . 4  Proposed  Mining  Method 

The  method  under  consideration  will  involve  multiple 
access  from  a  surface  mine  highwall,  and  headings  driven 
down  dip  (approximately  1500m  to  a  fault  boundary).  This 
concept  of  punch  mining  from  the  highwall  would  utilize  high 
capacity  equipment  to  extract  coal  from  the  portals  serving 
individual  mining  rooms.  As  the  parallel  rooms  are 
extracted,  and  the  portals  are  advanced  along  the  highwall, 
the  old  portals  become  redundant  and  will  be  covered  during 
reclamation . 


6 . 5  Stabi 1 i tv  Analysis 

A  number  of  possible  mine  geometries  have  been  studied 
with  the  assumed  rock  mass  parameters  and  initial  stress 
conditions  held  constant.  The  sensitivity  of  the  mine 
structure  to  variation  of  these  parameters  has  therefore  not 
been  considered.  The  stability  of  coal  pillars  and 
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surrounding  strata  have  been  considered  separately  for 
analytic  purposes,  but  interactions  at  roof  and  floor  will 
be  present.  The  geometries  analysed  represent  extraction 
ratios  of  60%,  55%  and  50%. 

6.5.1  Coa 1  Pillar  Strength 

The  cube  strength  of  coal  has  been  estimated  from  the 
uniaxial  compressive  strength  results.  Adjustments  have  been 
made  to  allow  for  sample  size  and  shape  using  the  following 
formula  based  on  Hustrulid  (1975): 

o  c  =  a  . -Je/ (0.778  +  0.222. W/H) 

ac  =  Compressive  strength  of  a  1  inch  cube 
act  =  Test  sample  compressive  strength 
H  =  Test  sample  height  (inches) 

W  =  Test  sample  width  (inches) 


Only  one  uniaxial  compressive  test  sample  result  was 
available  and  this  is  not  considered  to  be  sufficient  for 
normal  design  purposes,  but  with  no  other  data  available 
this  result  was  used.  For  the  test  sample  of  height 
2.02  inches  and  width  (diameter)  1.62  inches  the  uniaxial 
compressive  strength  was  15.6  MN/m2.  The  calculated  cube 
compressive  strength  is  21.0  MN/m2  for  a  1  inch  cube.  All 
mine  pillars  will  have  dimensions  greater  than  3  ft.  (0.9m) 
and  hence  the  pillar  cube  compressive  strength  will  be 
3.5  MN/m2  (one  sixth  of  the  one  inch  cube  strength).  The 
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size  factor  will  therefore  be  constant  and  the  shape  factor 
will  vary  according  to  the  following  equation  (Hustrulid, 
1975)  : 

a  =  a  (0.778  +  0.222  W/H) 
p  c 

Qp  =  Pillar  compressive  strength  (MN/m2 ) 

0Q  =  Cube  compressive  strength  (MN/m2) 

The  height  of  rooms,  and  hence  that  of  pillars  has  been 
assumed  to  be  2.5m  and  this  reduces  the  equation  to  the  form 
below: 


a  =  2.72  +  0.31  W 
P 

Cp  =  Pillar  compressive  strength  (MN/m2) 
W  =  Width  of  square  pillar  (m) 


The  width  used  in  the  equation  assumes  that  the  pillars  are 
square,  but  the  proposed  mining  method  will  not  utilise 
square  pillars.  Pillars  will  be  long  and  thin,  thus  an 
effective  width  should  be  used  for  design  purposes.  The 
effect  of  unequal  pillar  side  lengths  on  pillar  strength  has 
not  been  previously  studied  in  any  great  detail.  Sheorey  & 
Singh  (1974)  have  concluded  that  the  strength  of  a 
rectangular  pillar  is  the  same  as  that  of  a  square  pillar 
whose  sides  are  equal  to  the  average  of  the  two  rectangular 
pillar  widths.  This  result  indicates  that  for  a  given  pillar 
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area  a  rectangular  pillar  will  have  a  higher  strength  than  a 
square  one.  For  wide  pillars  with  confined  cores  this  may  be 
correct,  but  for  thin,  long  pillars  the  equivalent  square 
pillar  width  should  not  be  the  average  of  the  two  sides. 

This  is  understood  when  a  rectangular  pillar  with  one  side 
of  zero  length  is  considered.  The  equivalent  square  pillar 
will  have  a  width  equal  to  half  of  the  non-zero  length  of 
the  rectangular  pillar,  which  is  incorrect.  It  would  seem 
more  appropriate  to  base  the  design  of  a  rectangular  pillar 
on  its  equivalence  to  a  square  pillar  of  the  same  plan  area. 
Using  this  hypothesis,  and  the  assumption  of  a  pillar  length 
of  30m  parallel  to  the  long  axis  of  the  rooms,  the  pillar 
strength  relationship  becomes: 


W  =  Pillar  width  between  rooms  (m) 
a  =  Pillar  compressive  strength  (MN/mz ) 


This  relationship  is  presented  graphically  in  Figure  6.4.  It 
should  be  noted  that  this  relationship  is  only  valid  for 
pillars  with  a  height  of  2.5m  and  a  long  dimension  of  30m. 

6.5.2  Numer ica 1  Analysis 

The  DDSEAMS  program  was  used  to  analyse  three  layouts 
at  a  depth  of  150m  (Runs  21,  22  &  23)  with  the  graphical 
results  presented  in  Appendix  8.  The  mine  geometries  for 
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these  runs  are  shown  in  Figure  6.5  together  with  the  average 
vertical  stress  for  each  pillar.  The  rock  mass  and  coal  seam 
parameters  were  estimated  from  the  geotechnical  testing  and 
adjusted  for  j_n  s i tu  conditions  (tabulated  previously).  The 
overburden  rock  density  was  assumed  to  be  equivalent  to 
0.025  MN/m2  per  metre  of  depth.  No  data  for  the  pre-mining 
stress  field  was  available  so  hydrostatic  conditions  were 
assumed.  This  was  thought  to  be  reasonable  because  of  the 
shallow  nature  of  the  excavations. 

The  PSTRESS  program  was  run  on  the  data  from  the  three 
DDSEAMS  program  runs,  to  produce  values  for  the  stability  of 
the  surrounding  strata.  In  all  three  runs  the  roof  sandstone 
was  found  to  be  stable,  but  slight  instability  in  the  floor 
mudstone  was  observed.  The  stability  criterion  results, 

0.75m  into  the  floor  mudstone,  are  presented  in  Figure  6.7 
for  the  present  mining  room  and  the  previous  room. 

Pillar  safety  factors  (based  on  elastic  stress 
comparisons)  can  be  calculated  by  dividing  the  pillar 
strength  by  the  average  pillar  stress.  These  results  are 
presented  in  Figure  6.6,  where  pillar  strengths  of  6.12  and 
7.53  MN/m2  have  been  assumed  for  the  4m  and  6m  wide  pillars 
respect i ve 1 y . 


6 . 6  D i scuss ion  of  Resu 1 ts 

In  considering  the  stability  of  the  mine  structure  it 
will  be  assumed  that  only  the  room  being  mined,  and  the 
previous  room,  need  to  be  stable.  This  requires  that  the 
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abutment  pillar  (No.1),  and  pillars  No. 2  and  No. 3  should  be 
stable.  This  assumption  is  only  valid  when  failed  pillars  do 
not  exhibit  brittle  behaviour,  but  yield  progress i ve 1 y ,  thus 
avoiding  catastrophic  failures. 

6.6.1  Pillar  Stabi 1 i tv 

Using  the  tributary  area  theory,  the  expected  pillar 
stresses  for  extraction  ratios  of  60%,  55%  and  50%,  are  9.4, 
8.3  and  7.5  MN/m2  respectively.  For  these,  pillar-  widths  of 
15.4m  (60%  extraction),  10.8m  (55%  extraction)  and  7.8m  (50% 
extraction)  would  be  required  to  give  a  safety  factor  of  1. 

The  results  from  the  numerical  modelling  indicate  that 
the  proximity  of  the  abutment  to  the  last  two  pillars,  has 
the  effect  of  reducing  their  average  stress  values  below 
that  predicted  by  the  tributary  area  theory.  If  a  room  width 
of  6m  is  assumed,  then  the  tributary  area  theory  predicts 
that  a  pillar  width  of  7m  is  required  for  a  safety  factor  of 
1  (pillars  30m  long  and  extraction  ratio  of  46%).  The  pillar 
safety  factors  for  the  three  runs  (Figure  6.6),  predicted  by 
the  numerical  model,  indicate  that  extraction  ratios  above 
50%  should  be  possible.  The  results  illustrate  the 
undesirable  aspects  of  having  non-uniform  pillar  sizes,  due 
to  the  relatively  low  strengths  of  narrow  pillars.  As  the 
pillar  width  is  reduced,  the  strength  is  reduced,  and  the 
pillar  stress  is  increased.  This  indicates  that  maximum 
extraction  will  be  achieved  with  uniform  pillar  widths. 
Design  on  the  basis  of  a  safety  factor  of  1  is  considered  to 
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be  adequate,  as  the  worst  possible  case  of  greatest  depth  is 
ana  1 ysed . 

Pillar  stabilities  are  predominantly  controlled  by  the 
vertical  stresses,  which  can  be  reasonably  predicted  using 
overburden  density  assumptions.  The  effect  of  horizontal 
stress  upon  the  pillar  stability  is  not  expected  to  be 
great,  as  the  seam  is  gently  dipping.  The  assumption  that 
the  seam  can  be  analysed  as  being  flat  is  not  expected  to 
affect  the  results  by  more  than  5%,  and  the  approach  is 
therefore  considered  to  be  valid. 

6.6.2  Open i no  Stabi 1 i tv 

The  stability  of  the  roof  span  is  not  likely  to  be 
controlled  by  intact  rock  shear  failure,  as  it  consists  of  a 
massive  sandstone  strata  unit  over  6m  thick.  From  the 
limited  discontinuity  data  available,  steep  joint  sets 
(approximately  60°  dip)  are  known  to  exist,  but  their 
spacing  and  orientation  are  unknown.  Stability  of  wedges 
delineated  by  these  steep  joint  sets  and  bedding  planes  will 
be  greatly  affected  by  the  horizontal  stresses,  which  are  at 
present  unknown.  Consequently  no  complete  stability  analysis 
can  be  undertaken  at  present,  but  it  should  be  mentioned 
that  experience  has  shown  roof  beds  of  this  type  to  be  very 
stable  with  roof  bolting  support.  In  order  to  estimate  the 
support  requi rements ,  further  information  on  the  dip, 
orientation,  spacing  and  shear  strength  of  these 
discontinuities  is  required. 
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The  orientation  of  rooms  will  be  a  major  factor  in 
analysing  the  roof  and  rib  stability.  Slabbing  failure  of 
pillar  sides  into  rooms  is  dependent  upon  the  cleat 
direction,  spacing  and  continuity,  which  are  unknown.  Pillar 
slabbing  is  unsafe,  causes  production  delays,  and  also  has 
the  effect  of  increasing  opening  spans  and  reducing 
effective  pillar  strengths. 

Floor  stability  problems  are  often  encountered  when 
weak  floors  are  highly  stressed,  producing  failures  which 
migrate  into  the  entry  (floor  heave).  High  horizontal 
stresses  may  cause  buckling  failures  of  thin,  weak  floor 
beds,  which  are  further  pronounced  by  highly  stressed  pillar 
abutments  load-ing  the  floor  strata  beyond  its  bearing 
capacity.  The  effect  of  water  on  the  floor  may  be  a  major 
factor  for  the  mudstones  and  siltstones.  Swelling  Index  and 
Slake  Durability  tests  indicate  that  the  swelling  of  floor 
beds,  together  with  their  reduced  strength  when  wet,  will 
cause  problems  for  heavy  mobile  machinery.  Floor  heave  will 
probably  not  be  a  problem  as  it  is  a  time-dependent 
phenomenon  and  the  most  highly  stressed  zones  will  be  at  the 
ends  of  rooms.  The  stress  criterion  results  for  the  floor 
mudstones  (Figure  6.7),  based  on  geotechnical  testing  and 
stresses  from  the  numerical  modelling,  indicate  that  slight 
instability  of  floor  strata  at  the  centre  of  openings  will 
be  present . 

The  overall  analysis  has  been  based  upon  limited  test 
results  from  the  more  competent  portions  of  the  strata, 
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because  of  sample  preparation  problems  for  the  weaker  units. 
This  means  that  the  average  testing  results  will  be  biased 
towards  the  higher  strength  values.  Adjustments  have  been 
made  by  applying  scaling  factors  to  estimate  the  rock  mass 
parameters  based  on  previously  published  results.  These 
factors  are  not  well  understood  and  may  not  be  applicable  to 
the  strata  under  study,  thus  a  degree  of  uncertainty  has 
been  introduced. 

This  design  analysis  is  sufficient  for  the  initial 
technical  and  economic  feasibility  stages.  Better 
predictions  can  be  made  by  undertaking  jjn  situ  testing  of 
full  size  excavations,  at  sufficient  depth,  with  monitoring 
of  strata  behaviour. 


7.  CONCLUSIONS 


7 . 1  Model  1 inq  Techniques 

1 .  The  numerical  analysis  technique  adopted  provides 

improved  prediction  of  stress  distributions  compared  to 
simpler  techniques,  such  as  the  tributary  area  and  beam 
theory.  This  is  especially  true  in  the  proximity  of 
abutments,  for  workings  of  limited  lateral  extent  in 
relation  to  depth,  and  where  multiple  seam  interactions 
are  present . 


2 .  Layered  strata  with  stiffer  mechanical  properties  cannot 
be  modelled  by  assigning  unmined  coal  seams  in  the 
numerical  model  presented. 

3 .  Physical  modelling  techniques  provide  a  simple  means  of 
qualitatively  simulating  complex  failure  mechanisms  in 
mine  structures,  by  relating  their  behaviour  to  the 
failures  observed  in  real  mine  situations. 


4 .  The  use  of  a  stress  criterion  to  predict  zones  of 

overstressed  rock  is  useful  for  indicating  potential 
fail ure  zones . 


5 .  When  pillar  safety  factors  or  strata  stress  criterion 
values  predict  that  failures  will  take  place,  the 
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resulting  stress  re-distributions  are  not  accounted  for 
in  the  numerical  analysis. 


7 . 2  Multi -Seam  Studies 

1 .  The  numerical  model  predicted  zones  of  potential 
instability  that  correspond  to  observed  instability  in 
the  3-D  physical  model. 

2 .  The  pillar  failures  observed  in  the  3-D  physical  model 
were  similar  to  those  observed  in  the  underground 
phenomenological  study. 

3 .  The  base  friction  model  predicted  a  cave  angle  of 
approximately  70°  for  de-pi  1  Taring  operations  in  the 
lower  seam. 

4 .  The  numerical  model  only  predicts  potential  failure 
zones,  and  should  therefore  not  be  used  to  estimate  the 
stress  re-distributions,  as  it  does  not  model  these 
fail ures . 


7 . 3  Mine  Design  Simulation 


1 ,  From  geotechnical  testing  of  core  samples  of  the 
proposed  mine  strata  at  Coal  Valley  the  following 
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results  were  obtained: 

(a)  Uniaxial  Compression  Testing: 


Younq'  s  1 

Poisson' s 

Compress i ve 

Modulus 

Rat  io 

Strenath 

(MN/m2 ) 

(MN/m2 ) 

SANDSTONE  1.36  x  IQ4 

0.267 

73.9 

(3  samples) 

SILTSTONE  1 .00  x  104 

0.288 

66.2 

(2  samples) 

MUDSTONE  8.62  x  103 

0.172 

70.2 

(1  sample) 

COAL  *  1 .92  x  103 

- 

15.6 

(1  sample) 

*  undercored  sample  (no  strain 

gauges  attached) 

(b)  Brazilian  Testing  (tensil 

e  strength): 

SANDSTONE  (14  samples) 

=2.29  MN/m 

2 

SILTSTONE  (6  samples) 

=  3.47  MN/m2 

MUDSTONE  (11  samples) 

=  3.97  MN/m2 

COAL  (4  samples) 

=  0.91  MN/m2 

(c)  Triaxial  Compression  Testing: 

(confining  stresses  less  than  15  MN/m2) 
Linear  Mohr /Cou lomb  failure  envelopes: 


SANDSTONE:  t  =  6.9  +  crTan4  5 

MUDSTONE:  t  =  6.9  +  aTan38° 

COAL:  T  =  3.0  +  aTan45° 

t  =  Shear  stress  (MN/m2) 
a  =  Normal  stress  (MN/m2) 

(d)  Direct  Shear  Testing  of  Discontinuities: 

MUDSTONE  JOINT  (50°  dip):  Friction  angle  =  19.9° 
MUDSTONE  JOINT  (30°  dip):  Friction  angle  =  27.1° 
MUDSTONE  BEDDING  (15°  dip):  Friction  angle  =  23.3 
COAL  BEDDING  (15°  dip):  Friction  angle  =  23.8° 

(e)  Slake  Durability  Testing  of  Footwall  Mudstones 

SDI  (average  of  3  samples)  =  63.5% 

SDI  =  Slake  Durability  Index  (second  cycle) 


(f)  Swelling  Index  Testing: 

SSI  (average  of  2  samples)  =  2.7% 

SSI  =  Swelling  Strain  Index  (unconfined) 


' 
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2 .  The  estimated  rock  mass  parameters  from  the  empirically 
adjusted  geotechnical  results  are: 


Mater ia 1 

E  (MN/m2) 

G  (MN/m2) 

ROE 

ROCK 

3.9  x  103 

1.6  x  103 

0.19 

COAL 

3.9  x  102 

1.5  x  102 

0.30 

E  =  Young's  Modulus 
G  =  Shear  Modulus 
ROE  =  Poisson' s  Rat io 


3 .  The  estimated 
the  rock  mass 

Material 

SANDSTONE 

MUDSTONE 

COAL 


linear  Mohr /Coulomb 
types  are: 

C  (m/m2) 

2.3 

2.3 

0.55 


failure  envelopes  for 


PHI  (degrees ) 

40.0 

33.0 

40.0 


C  =  Rock  Mass  Cohesion 

PHI  =  Internal  Friction  Angle 


4 .  The  estimated  pillar  strength  for  long  pillars  (30m 
parallel  to  rooms)  is  given  by: 
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c  =  2.72  +  1 . 7/w~ 

P  t 

W  =  Pillar  width  between  rooms  (m) 

Qp  =  Pillar  compressive  strength  (MN/m2 ) 

5 .  This  initial  analysis  indicates  that  extraction  ratios 
(by  area)  greater  than  50%  may  be  possible  at  150m 
depth . 

6 .  The  analysis  is  based  on  geotechnical  testing  of  core 
from  shallow  seam  conditions  which  may  not  be 
representative  of  the  strata  conditions  at  greater 
depth . 


8.  RECOMMENDATIONS 


8 . 1  Model  1 inq  Techniques 

1 •  Predictions  of  actual  rock  mass  behaviour  and  failure 
mechanisms  should  be  improved.  This  could  be  done  by 
comprehensive  geotechnical  testing  of  intact  rock  and 
discontinuities,  definition  of  structural  geology 
(including  joint  surveys)  and  measurement  of  j_n  situ 
stresses.  Monitoring  of  the  defined  rock  structure 
should  be  carried  out  during  excavation,  and  back 
analysis  techniques  used  to  provide  an  accurate  means  of 
predicting  the  behaviour  of  the  mine  structure  from  data 
collected  prior  to  mining. 

2 .  Numerical  methods  should  be  developed  to  simulate  the 
actual  rock  mass  behaviour  when  failures  occur.  This 
could  possibly  be  achieved  by  combinations  of  boundary 
element,  finite  element  and  block  models.  An  example  of 
this  could  be;  to  define  an  active  caving  zone  (caved 
rock)  by  means  of  a  block  model,  this  could  interface 
with  a  finite  element  model  from  the  active  caving  zone 
to  the  caving  limits  (failed  rock),  which  is  then 
interfaced  with  a  boundary  element  model  for  the 
remainder  of  the  rock  mass  (pre- f ai lure ) . 

3 .  Pillar  failure  mechanisms  should  be  investigated 
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further,  under  actual  mine  conditions,  to  improve  pillar 
strength  prediction  capabilities.  This  should  include 
time  dependent  analysis  as  pillar  deterioration  with 
time  is  commonly  observed  in  coal  mines. 

4 .  Where  significant  failures  are  present  in  the  mining 
structure,  elastic  analysis  techniques  do  not  provide 
adequate  estimation  of  the  resulting  stress 

re-di str ibut ions .  Therefore  it  is  recommended  that  more 
complex  analysis  techniques,  that  model  these  failure 
mechanisms,  should  be  employed. 

5 .  The  limitations  of  modelling  techniques  should  be 
investigated  further,  as  the  ability  to  model  actual 
rock  mass  behaviour  is  still  in  its  rudimentary  stage. 
Better  techniques  are  constantly  being  developed,  but 
their  applicability  in  modelling  a  complex  reality 
contain  many  uncertainties. 

6 .  The  applicability  of  two-dimensional  analysis  techniques 
should  be  investigated  further,  as  they  are  often  used 
because  three-dimensional  techniques  are  too  complex  and 
expensive  to  undertake. 
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8 . 2  Multi -Seam  Studies 


1 .  Where  caving  phenomena  are  encountered,  analysis 
techniques  that  model  these  complex  failures  should  be 
emp 1 oyed . 

2 .  The  rock  mass  behaviour  and  multi -seam  interactions 
during  caving  should  be  investigated  further  to  enable 
realistic  modelling  of  mine  structures. 

3 .  Coal  pillar  strength  relationships  should  be 
investigated  further  to  improve  prediction  capabilities 
for  Rocky  Mountain  conditions. 


8 . 3  Mine  Design  Simulation 

1 .  Further  geotechnical  testing  of  strata  conditions  at 
greater  depth  should  be  carried  out. 

2 .  Monitoring  of  the  rock  mass  behaviour  around  actual  mine 
openings,  at  the  depth  conditions  expected,  should  be 
undertaken . 


3 .  Time-dependent  behaviour  of  mine  pillars  and  strata 
should  be  investigated,  as  the  proposed  mining  method 
will  probably  involve  pillars  that  exhibit 
time-dependent  yield  behaviour. 
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More  complex  modelling  techniques  should  be  employed  for 
prediction  of  the  behaviour  of  mine  structures.  The 
results  of  strata  monitoring  for  experimental  openings 
would  be  required  for  realistic  modelling. 
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PILLARS 


FIGURE  2.1  SQUARE  PILLAR  LAYOUT  PLAN 
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FIGURE  3.1  CONFIGURATION  OF  FLAT  ELLIPTICAL  CRACK 


VERTICAL  STRESS  (lb/ft 
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FIGURE  3.2  ANALYTICAL/NUMERICAL  VERTICAL  STRESS  COMPARISON 
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FIGURE  3.3  SIMPLIFIED  MINE  LAYOUT 


RUN  NUMBER 
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FIGURE  3.4  TABLE  OF  MAJOR  PARAMETERS  FOR  RUNS  06  TO  11 


SHEAR  STRESS  (MN/mV  HORIZONTAL  STRESS  ( MN/ m‘ )  VERTICAL  STRESS  (  MN/mO 
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VERTICAL  STRES S  IN  S  E AM 


SHEAR  STRESS  (1.5m  INTO  ROOF) 


FIGURE  3.5 


STRESS  RESULTS  FOR  RUN  08 


RUN  NUMBER 
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FIGURE  3.6  TABULATED  COMPARISON  OF  SIMPLIFIED  LAYOUT 


86 


Run  14 


DEPTH 
1000  m 


Run  1 5 


FIGURE  3.7  PI LEAR  VERTICAL  STRESSES 


RUNS  14  &  15 
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Run  1 7 
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FIGURE  3.8  PILLAR  VERTICAL  STRESSES 


RUNS  17  &  18 


SHEAR  STRESS 
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FIGURE  4.1  MOHR/COULOMB  FAILURE  CRITERION 
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Figure  5.1  Mine  Layout  for  Run  19 
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Figure  5.2  Numerical  Results  for  Upper  Seam  (Run  19) 
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Figure  5 . 3  Numerical  Results  for  Lower  Seam  (Run  19) 


Vertical  Stress  (  MN/m2 ) 
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(w)  ludea 


Horizontal  Distance  (m) 


Horizontal  Stress  ( MN/m2) 
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(UJ)  ludea 


Figure  5.5  (Run  19)  Horizontal  Distance  (m) 


Major  Principal  Stress  (  MN/m2 ) 
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(m)  Miclsa 


Figure  5.6  (Run  19)  Horizontal  Distance  (m) 


Minor  Principal  Stress  (  MN/m2 ) 
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Figure  5.7  (Run  19)  Horizontal  Distance  (m) 


Maximum  Shear  Stress  (  MN/m2 ) 
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(tu)  ludea 


Figure  5.8  (Run  19)  Horizontal  Distance  (m) 


STRESS  CRITERION  (Negative  Areas  Shaded) 
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FIGURE  5.10  BASE  FRICTION  MODEL  LAYOUT 
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Lower  seam 


Upper  seam 


FIGURE  5.1 


3-D  PHYSICAL  MODEL  LAYOUT  PLAN 
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3-D  PHYSICAL  MODEL  (RUN  1) 


NUMERICAL  MODEL  (RUN  19) 


FIGURE  5.12 


NUMERICAL 


/  3-D  PHYSICAL  MODEL  COMPARISON 
(OPENING  2) 
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SIMPLIFIED 
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FIGURE  6.1  TRIAXIAL  TEST  RESULTS  -  SANDSTONE 
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FIGURE  6.2  TRIAXIAL  TEST  RESULTS  -  MUDSTONE 
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FIGURE  6.4  PILLAR  STRENGTH  AS  A  FUNCTION  OF  WIDTH 


"All  seam  geometries  at  150m  depth  and  2.5  m  thickness 
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FIGURE  6.5  MINE  SIMULATION  LAYOUTS  AND  PILLAR  VERTICAL  STRESSES 
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FIGURE  6.6  MINE  SIMULATION  LAYOUTS  AND  PILLAR  SAFETY  FACTORS 
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FIGURE  6.7  STRESS  CRITERION 


(0.75m  INTO  FLOOR) 
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Plate  1.  BASE  FRICTION  APPARATUS 


Plate  2.  MTS  SERVO  CONTROLLED  STIFF  TESTING  MACHINE 
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Plate  3.  3-D  MODEL  IN  COMPRESSION  FRAME 
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Plate  4.  BASE  FRICTION  MODEL  -  CAVING  INITIATED 
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Plate  5.  BASE  FRICTION  MODEL  -FINAL  CAVE 
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Plate  6.  3-D  MODEL  (RUN  1)  -  YIELD  ST£TE 
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Plate  7.  3-D  MODEL  (RUN  1,  OPENING  No . 2 )  -  YIELD  STATE 
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Plate  8.  3-D  MODEL  -  AFTER  TEST 


(RUN  1,  PLAN  VIEW,  UPPER  SEAM) 
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Plate  9.  3-D  MODEL  (RUN  2,  LOWER  SEAM)  -  BEFORE  TEST 


Plate  10.  3-D  MODEL  (RUN  2,  UPPER  SEAM) 


BEFORE  TEST 
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Plate  11.  3-D  MODEL  (RUN  2,  SIDE  VIEW)  -  BEFORE  TEST 


117 


Plate  12.  3-D  MODEL  (RUN  2,  SIDE  VIEW)  -  0.78  MN/m 
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Plate  13.  3-D  MODEL  (RUN  2,  SIDE  VIEW)  -  YIELD  STATE 
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Plate  14. 


3-D  MODEL  -  AFTER  TEST 

(RUN  2,  PLAN  VIEW,  LOWER  SEAM) 
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Plate  15.  3-D  MODEL  (RUN  2,  PLAN  VIEW,  UPPER  SEAM)  -  AFTER  TEST 
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Plate  16  *  NEWLY  EXPOSED  PILLAR  SIDE  AND  SUPPORTS 
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Plate  17.  PILLAR  EXPANSION  NEAR  ROOF 
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Plate  18.  PILLAR  FRACTURING 
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Plate  19.  PILLAR  EXPANSION  AND  FAILURE  OF  SUPPORTS 
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Plate  20.  PILLAR  EXPANSION ,  FLOOR  HEAVE  AND  SUPPORT  FAILURE 
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Plate  21.  CAVED  AREA  AND  FLOOR  HEAVE 
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Plate  22.  CAVE  OF  IMMEDIATE  ROOF  BEDS 
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Chapter  2 
INTRODUCTION 

The  design  of  underground  mining  structures  for  tabular 
orebodies,  such  as  coal  seams,  is  primarily  controlled  by 
stability  analyses.  In  order  to  make  an  assessment  of  the 
stability  of  an  excavation  it  is  necessary  to  estimate  the  total 
stresses  and  displacements  after  excavation  has  taken  place  (the 
total  stresses  are  the  sum  of  the  initial  stresses  and  those 
induced  by  the  excavation,  whereas  the  displacements  are  only 
those  induced  by  the  excavation  i.e.  the  initial  displacements 
are  zero).  These  stresses  and  displacements  are  then  compared 
with  a  set  of  stability  or  failure  criteria  to  assess  the 
stability  of  the  proposed  design. 

The  total  stresses  and  displacements  after  excavation  has 
taken  place  are  functions  of:  (i)  the  initial  stress  state  prior 
to  excavation,  (ii)  the  geometry  of  the  excavation,  (iii)  the 
rock  mass  characteristics,  and  (iv)  the  excavation  sequence. 

(i)  The  initial  stress  state  is  usually  unknown  and  must  be 
measured  or  estimated  from  knowledge  of  similar  strata  conditions 
e 1 sewhere . 

(ii)  The  excavation  geometry  is  a  three  dimensional  surface  which 
is  usually  complex  in  reality  and  thus  is  normally  simplified  to 
permit  analysis  to  be  undertaken  on  a  reasonable  scale. 

(iii)  The  rock  mass  characteristics  define  the  relation  between 
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the  stress  changes  and  the  displacements  that  they  induce,  the 
nearer  this  relation  is  to  reality,  the  greater  the  complexity  of 
analysis  required  for  solution.  The  relation  may  be  elastic 
(linear  or  non-linear),  plastic  (rigid/ideally  plastic, 
e 1  as t i c/ i dea 1 1 y  plastic  or  elasto-plastic) ,  or  creep 
(visco-elastic  or  visco-plastic).  The  material  properties  may  be 
the  same  in  all  directions  (isotropic),  equal  in  two  directions 
(transversely  isotropic)  or  different  in  all  three  directions 
(orthotropic).  Geologically  the  rocK  strata  usually  has 
discontinuities  (joints  and  faults)  and  the  analysis  may  take 
this  non -homogeneous  rock  mass  into  account. 

(iv)  The  excavation  sequence  is  only  important  where  the  rock 
mass  characteristics  have  irrecoverable  displacements  as  in 
plasticity  (yield  functions),  creep  (time  dependent  yield)  and 
slip  surfaces  (joints,  faults).  These  relationships  exhibit  path 
dependent  effect  upon  the  final  results  and  require  that  the 
excavation  sequence  be  examined  as  well. 

Historically  the  types  of  analyses  that  were  possible,  were 
simple  and  very  restrictive  in  their  applications.  Usually  the 
geological  structure  had  to  be  ignored,  and  the  excavation 
geometry  and  material  behaviour  greatly  simplified  in  order  to 
provide  closed  form  solutions.  With  the  advent  of  computers, 
numerical  methods  have  become  available  which  enable  complex 
geometries  in  three  dimensions,  with  complex  material  behaviour 
and  the  modelling  of  geological  reality  to  be  undertaken. 
However,  these  complex  formulations  are  both  difficult  to  set  up, 
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and  expensive  to  run,  which  usually  precludes  them  from  use  in 
practical  mine  design. 

The  errors  involved  in  quantifying  the  rock  mass 
characteristics  by  measurement  or  estimation  from  intact  material 
properties  and  geological  structure  are  significant  enough  to 
question  the  use  of  complex  rock  behaviour  assumptions  in  the 
analysis.  It  has  been  shown  that  for  relatively  deep  excavations 
the  zone  of  inelastic  behaviour  (region  of  yielding  and  failure) 
is  usually  only  confined  to  the  immediate  vicinity  of  the 
excavation  and  that  the  remainder  of  the  rock  mass  may  be  assumed 
to  behave  elastically. 

There  exists  the  need  for  a  method  of  analysis  of  mine 
excavations  which  is  based  on  relatively  simple  assumptions  of 
rock  mass  characteristics,  is  able  to  simulate  sufficiently 
complex  excavation  geometries,  and  is  both  simple  and  cheap  to 
run,  thus  enabling  many  layouts  to  be  analysed.  This  usually 
requires  that  plane  strain  formulations  are  used  which  reduce  the 
problem  to  a  two  dimensional  one,  and  the  assumption  of  linear 
elastic  behaviour  of  the  rock  mass  to  simplify  the  analysis. 
These  simplifications  are  adequate  when  the  main  requirement  is 
for  an  analysis  technique  which  is  used  to  examine  the 
sensitivity  of  the  overall  design  with  respect  to  variation  of 
each  of  the  major  parameters  in  turn. 

Finite  element  and  finite  difference  computer  programs 
divide  the  whole  area  of  influence  of  the  excavation  into 
discrete  elements  and  approximations  are  made  to  each  element. 
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This  makes  them  relatively  powerful  for  modelling  complex 
geometries  and  material  properties,  but  time  consuming  to  prepare 
and  costly  to  run.  Boundary  element  programs  only  divide  the 
boundary  surfaces  into  discrete  elements,  from  which  the  stresses 
and  di spl cements  at  any  point  in  the  rock  mass  can  be  calculated. 
Because  of  the  reduced  number  of  elements  with  boundary  element 
programs  in  comparison  to  finite  element  programs,  they  are 
usually  much  easier  to  set  up  and  cheaper  to  run. 

The  DDSEAMS  program  described  in  this  report  is  a  boundary 
element  program  which  utilises  the  Displacement  Discontinuity 
Method  as  developed  by  S.L. Crouch.  This  program  is  a 
two-dimensional  formulation  (plane  strain)  for  analysis  of 
multiple  seam  extraction  of  linearly  elastic  coal  seams  in  an 
orthotropic  linear  elastic  half-space.  The  seams  are  assumed  to 
have  negligible  width  in  comparison  to  their  length  and  lie 
parallel  and  not  too  close  to  the  surface.  The  program,  as 
developed  within  the  department,  can  handle  up  to  five  seams  with 
up  to  eighty  elements  per  seam,  and  is  suitable  for  both 
interactive  and  batch  operating  modes.  A  special  feature  of  the 
program  is  that  it  can  produce  graphical  output  of  stresses  and 
displacements  within  the  area  of  interest  whilst  in  interactive 
mode . 
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Chapter  3 

MATHEMATICAI  THEORY 

This  section  will  only  summarise  the  major  stages  in  the 
solution  process  and  all  detailed  equations  and  some  intermediate 
steps  are  shown  in  re f.  1.  The  numerical  solution  procedure 
is  described  in  the  next  section,  and  this  section  will  only  deal 
with  the  derivation  of  the  influence  functions  for  the 
fundamental  displacement  discontinuity  line  segment  that  is  used 
in  the  solution. 

The  theory  of  elasticity  is  used  to  calculate  the  effect  of 
a  constant  displacement  discontinuity  over  a  finite  line  segment 
parallel  to  the  surface  of  a  semi -inf  ini te,  homogeneous, 
orthotropic  linear  elastic  body  y>0  (see  Fig.  3.1). 
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The  special  case  in  which  the  directions  of  elastic  symmetry  of 
the  material  are  parallel  to  the  x  and  y  axes  will  be  analysed, 
and  plane  strain  conditions  will  reduce  the  problem  to  a 
two-dimensional  one.  The  following  boundary  conditions  are 
assumed  : 

(i)  The  displacements  are  continuous  everywhere  in  the  body, 
except  over  the  line  segment  in  question. 

(ii)  The  normal  and  shear  stresses  are  zero  all  along  the  surface 
of  the  half-space  (ground  surface). 

(iii)  All  stresses  and  displacements  are  zero  at  infinity. 

The  problem  considered  will  be  the  solution  for  the  line 
segment  shown  in  Fig.  3.1  where  : 


|x|  <  a  ,  y  -  -H  ( h  >  o ) 


These  conditions  can  be  expressed  as 


Where 


O'*;,  =  <  X  <  CO  ;  y  -  Q  , 

The  subscripts  on  (h.andlij  denote  the  negative  and  positive  sides 


of  the  line  y  =  -h  (i.e.  h^denotes  the  top  of  the  crack) 

and  D^are  the  relative  shear  and  normal  displacements  between 
the  two  surfaces  of  the  crack. 
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Uxand  are  the  actual  displacements  of  the  cracK  surfaces  in  the 
x  and  y  directions  respectively. 

The  solution  to  this  problem  is  found  using  the  method  of 
images  (Ref.  1),  giving  the  results  in  the  form: 


A  I  S 


Ut  = 

Ui  +  UL  + 

u. 

11 

:p 

b 

A  I 

cr.  +  cr.  + 

cj  uj 

S 

cr.. 

*0 

ii 

.  j_ 
'Ll? 

Where  the  superscripts  denote  the  solution  types: 

A  =  Actual  displacement  discontinuity 

I  =  Image  displacement  discontinuity 

S  =  Supplemental  displacement  discontinuity  solution 

The  index  notation  refers  to  the  usual  tensorial  notation,  where 

the  direction  of  the  x  and  y  axes  are  replaced  by  xi and  xl 

respectively,  thereafter  the  subscripts  1  and  2  refer  to  the 

direction  in  which  the  function  they  define  acts. 

The  solution  is  obtained  from  the  solution  of  a  constant 
displacement  discontinuity  over  a  finite  line  segment  in  an 
infinite  body  and  superimposing  an  image  solution  to  create  a 
shear  stress  free  surface  (ground  surface),  and  a  supplemental 
solution  to  reduce  the  normal  stresses  to  zero  along  the  same 
surface . 
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3.1  Actual  Displacement  Discontinuity 


The  solution  for  the  actual  constant  displacement 
discontinuity  over  the  line  segment  in  an  infinite  body 

Dl  —  C  Dae  }  Dy)  ^  |x|<CL  }  Lj  -  -  h 

A  A 

can  be  written  as  follows  for  U*  and  cr^:— 


-a 


_ L.^£k>»A) ! 

0+2.1  ^3,  (»*f  J  ^3*. 


D, 


*0 


JL  -  (*.,  j.-A.) 

~1>  1  *X  (l-fx) 


A 


a 


xx 


0,0. 


J.  _  1  Yf(*.y,*Q 


y**0 


Where 


-  O  +  g.Ki  +  gi) 

A"  (2,-2l) 


qrcban^v+K^  _  fy;4h:)  flrcbn 


—  ^3c'al^i«/{x-a)  +  (x *a) ^j0c*af+T^7hif 


y  =  scaled  y  co-ordinates 

0  =  Vir, 

}f-v=  solution  constant 

?i,=  Cj/Cc^  +  Cj  (t=.,x) 


C ,,  ,  Cu ,  C(ietc.  are  the  material  constants  for  relating  stresses 
to  strains. 
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The  full  equations  can  be  seen  in  ref.  1 


3.2  Image  Displacement  Discontinuity 


i  r 

Expressions  for  the  displacements  U^and  stresses  qjdue  to 
the  image  displacement  discontinuity  in  an  infinite  body  (shown 

*  r 

in  Fig.  3.2)  may  be  written  in  a  similar  way.  e.g.  for  U^and  cr^:- 


c|0+i.)  iy, 


D 


ypt.) 


Ax 

*,yi->0 


Where 


_  0 

Z-rr(t,-: 


—  (x-o-)  J(x-e.)\(y..^y  +  (x  +  a) 
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3.3  Supplemental  Displacement  Discontinuity 

By  superposition  of  the  actual  and  image  stresses  on  the 
surface  y  =  0  ,  the  normal  stresses  are  : 


flfi 

A 

r 

X 

ii 

KX 

+ 

and 

the  shear 

stresses  are 

zero  : 

R 

X 

A 

I 

+  CTxy  = 

°*3  = 

- 

The 

norma  1 

stresses  are  eliminated  by 

superimposing  a 

supplemental 

di spl acement 

discontinuity  for  the 

half  space  y<0 , 

that  has  appropriate  stress  boundary  values  on  y  =  0  : 


S 


Oku 

=:  O 

D 

5 

X 

-  <  :>c  <  +  tj  = 

O 

o-y3 

= 

The  boundary 

s 

value  of  cr33 

is  defined  in  terms  of 

the 

image 

di splacement 

di scont i nui ty 

because  this  solution 

does 

not 

introduce  additional  discontinuities  in  the  half  space  y^O . 


The  supplemental  solution  may  be  written  in  terms  of  a 
single  potent i a  I  ,X(  x , y ] 


i  a  1  ,^E(  x  ,  y )  and  the  stresses  ax„and  cr„are  found  : 


*3 


=  -c 


[Vi  .  Vjlo.a*)] 

64  A^bi  J 

Py,  b7  3?  ~| 

a w  -  "c“|_  J 
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Where  : 


2  S 

V  3E(*,y) 


b) 


o 


so  that  cryis  automatically  zero  when  y  =  0. 

6  s 

This  can  be  reduced  to  the  following  for  U  and  cr 

=7  Xac 


u  =- 

x. 


Z  p, 

VX 


~h  ^  ^ 

VX 


1  AiC*,b,-0  _  ^£6t,arKi')l  _  X  fW*,a,-iO  _ 

^y,  j  ay*  J 

2L a.-X)-  X i£(Xy,-*vi)l  -  X  ,) _  X ‘HfC*, 

l  +  ?l\  ^3C  S  3C  J  UJjL  ^3C  J 


Crx 


2C(iDx 


*,-X 


2C«Db 


i  f^Vfx3a,-h^_^V('jc,  _  jjVfC 

****.  J  Xl***; 


cxj^-hoT. 

iu  **  ^y,.  J 


X  -X 


^  •f  C*i Si~T»^ _  x ^  -f  _ L  3X  ^  _  X  ^  ~f  fe,  I 


X  L  X,1 


"a,1 


M  >a* 


>bx* 


Where  f(x,y^  -hL  )  has  the  same  meaning  as  for  the  image 
displacement  discontinuity  solution. 


3.4  Complete  Solution 

The  complete  solution  to  the  problem  is  given  by  the  sum  of 
the  three  separate  components,  using  the  principle  of 
superposition.  (The  full  equations  can  be  seen  in  ref.  1 .  ) 
The  equations  are  elaborated  by  substitution  of  the  appropriate 
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derivatives  of  the  function  f(x,y).  The  results  can  be  presented 
in  the  following  manner  : 


O'**  —  y)  ~  G  (x-a.)  +  ^G-fx-w^y)  -  G 


(x+a,  y)  _  G  (x-cy  y)JlG  +  £ G  (x+<x,  ~  £  Cx-a,  y)j  2>y 

$  ~  ^X+A,y)  -  G  ('x-a^y'jj  j)*  -f  G  (Tx  +  a(y)  -  G  fx-a 


u*  = 


—  j^G  (x*a,  <j)  -  G  (x-at  y  )j  1)^.  4-  [<  -  G  (x-a^y)j  7)t 


Ulj  =  j^G  -  G  C 


x-a 


a) 


Dx  + 


G  (x+atLy)  —  G  Gc-a^y)  T)t 


I  o 


lo 
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Chapter  4 

NUMERICAL  PROCEDURE 

The  mathematical  theory  in  section  3.  has  only  discussed  a 
constant  displacement  discontinuity  over  a  finity  line  segment. 
The  distribution  of  stresses  and  displacements  within  the 
structure  has  been  calculated  as  a  function  of  the  closure  (roof 
to  floor  convergence)  and  ride  (shear  displacement  between  roof 
and  floor)  components.  In  reality  the  closure  and  ride  components 
of  an  excavation  within  a  seam  are  not  constant  over  the  whole 
span,  but  are  functions  of  the  position  along  the  opening. 


4.1  Single  Seam 

This  variation  of  displacement  discontinuities  can  be 
numerically  modelled  by  treating  the  excavation  as  a  series  of 
finite  line  segments  joined  end  to  end  as  in  Fig.  4.1 
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Using  the  principle  of  superposition  the  total  induced 
stresses  and  displacements  at  any  point  are  found  by  sunmation  of 
the  effect  of  each  displacement  discontinuity  element  (seam 
element).  The  only  problems  that  remain  are  the  treatment  of 
mined  and  unmined  seam  elements,  the  calculation  of  the  closure 
and  ride  components  (D  and  D  )  to  be  applied  to  each  seam 
element,  and  the  treatment  when  complete  closure  between  roof  and 
floor  takes  place.  The  seam  is  defined  by  assuming  that  each 
element  which  has  not  been  mined  has  a  normal  stiffness  K  and  a 

t\ 

shear  stiffness  K  . 


The  induced  stresses  at  the  centre  of  the  i  th.  seam  element 
may  be  expressed  in  terms  of  the  displacement  discontinuities 

l  «r 

(closure  D  and  ride  D  ,  i  =  1 , N )  at  all  N  elements  as  follows  • 

ft  s 


c 

Cr 

s 


L 

a 

n 


where  fl  ,  fl  ,  R  R 


(i=  >,n) 


EQ.  1+-.1 


5»  JW  ns 

are  the  influence  coefficients  relating  induced  stresses  due  to 


each  seam  element 
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R 


X  -  X  +  (X  } 


x  -  x.  +  a 


L  J 


'  L  J 

x  -x  +  a. 


x  +-  a 


j 


-K)  - 


-K)  _  G  (  x  -  5c  -a  -J,  ) 

It-  1 


h)  -  G  ( 
6  V 


X  -  x  -  a  , 


-O 


-t.) 


£Q.  Jf.2. 


where  the  G  functions  are  obtained  from  ref  1. 

It  is  assumed  that  the  initial  stress  state  if  such  that  the 
principal  stresses  are  in  the  direction  of  the  x  and  y  axes  and 
that  the  magnitude  of  the  vertical  stress  is  proportional  to  the 
rock  density  and  the  depth. 


i  .  e 


A  system  of  2N  simultaneous  linear  equations  is  formed  for 
the  required  solution  of  2N  displacement  discontinuities  acting 
on  N  seam  elements  by  equating  EQ  4.1  and  EQ  4.3.  This  solution 
will  only  be  valid  if  all  the  elements  are  mined,  and  there  is  no 
total  closure  between  roof  and  floor. 

To  simulate  an  unmined  element,  that  element  is  assigned 
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induced  stress  values  which  are  proportional  to  the  displacement 
discontinuity  values  : 


Cr  _  K  D 
s  s  s 


6-  =  k  d 

k  n  n 


(  c  =  I,n) 


E  Q.  4-.1+- 


"The  final  equations  are  : 


EQ.  4-5 

MINED 

UNMI  N  ED 

hi  M 

)  A  D  +  )  A  D  ~ 

l_  SS  s  [_  }n  * 

j  =  1  • 

O 

K  D 

s  s 

N  N 

r  ^  -  v  ‘j  j 

)r  d  +  )  a  d  — 

/_V\  Si  i_  nn  n 

j - 1  j  =  i 

—  |l 

o 

k  D 

n  v\ 

where  each  element  is  defined  as  mined  or  unmined  and  the 
relevant  equation  is  associated  with  it. 


4.1.1  Complete  Closure 

So  far  no  account  has  been  taKen  of  the  possibility  of 
complete  closure  taking  place  (a  common  phenomenon  in  coal  mining 
methods)  and  this  situation  is  accounted  for  during  the  iterative 
solution  of  the  series  of  equations  in  EQ  4.5  by  the  following 
steps  : 

(1)  At  each  stage  of  iteration,  the  displacement  of  seam  elements 
approaches  more  closely  the  final  solution  (the  iterative 
procedure  is  described  later).  These  values  may  be  denoted  as 

i  i 

D(k)  and  D  (k)  for  the  i  th.  element  after  the  k  th.  iteration 

"  5 
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cycle. 

(2)  If  the  i  th.  element  is  mined,  the  closure  D  (K)  is  tested  to 
see  whether  it  is  greater  than  or  equal  to  the  seam  thickness.  If 
this  is  the  case,  we  assume  that  the  solution  for  this  element  is 
given  by  : 

D  =  o  ,  D  r  h0 

S  n- 

i.e.  no  shear  stress  is  transferred  between  the  roof  and  floor 
and  the  closure  is  held  equal  to  the  seam  thickness  h  . 

O 

(3)  After  the  next  cycle  of  iteration,  the  stresses  acting  across 
all  elements  that  were  assumed  to  undergo  complete  closure  during 
the  previous  cycle  are  computed.  If  the  normal  stress  (cr  )  is 

n 

tensile  (negative)  at  any  such  element,  then  the  complete  closure 
restraint  is  removed  and  the  normal  displacement  discontinuity  is 
allowed  to  change. 

Because  the  roof  and  floor  are  allowed  to  slide  freely  over 
each  other  (zero  shear  stress)  then  the  solution  is  path 
independent  (i.e.  independent  of  the  mining  sequence)  and  a 
unique  solution  is  obtained.  Once  the  solution  is  found,  the 
stresses  (hence  strains)  and  displacements  at  any  point  may  be 
found  by  calculting  the  influence  coefficients  of  each  element 
and  using  these  combined  with  the  displacement  discontinuity 
solut ions . 
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4.1.2  Iteration  Process 


The  2N  simultaneous  linear  equations  can  be  written  in  the 


form  : 


where  (no  closure) 


m  l  N  E  D 

UNin  /  NJ  E  D 

L 

L- 

B  - 

o 

k  D 

s 

s  s 

L 

B 

—  U  Yr 

k  D 

Kl 

n  * 

If  the  i  th.  equation  has  experienced  complete  closure,  then 
the  i  th.  seam  element  conditions  are  defined  and  the  i  th. 
equations  are  eliminated  from  the  system  provided  the  next  cycle 
does  not  provide  tensile  stresses  in  that  element,  thus  causing 
it  to  be  freed  again. 


The  system  of  equations  is  solved  by  Gauss-Seidel  iteration 
with  over -rel axa t i on .  The  Gauss-Seidel  iterative  procedure  is  a 
method  where  an  assumption  is  first  made  for  the  solution  values 

L  ^ 

of  D  and  D  and  during  each  iteration  cycle  a  better  approximation 
n.  s 
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for  the  value  of  each  unknown  is  calculated.  The  method  is  made 

more  efficient  by  using  values  which  have  been  calculated  in  the 

present  cycle  for  the  calculation  of  the  remainder  of  unknowns  in 
that  cycle.  A  further  optimisation  is  achieved  by  using 

over-relaxation,  where  convergence  of  the  solution  is  accelerated 
by  calculating  the  convergence  (difference  in  solution  value 
since  the  previous  cycle)  of  each  unknown  and  applying  a 
weighting  factor  (over-relaxation  factor)  to  give  a  better 

approximation  : 

=  (D)'"\  <j  (A  ?)'""> 

where  : 

GJ  ~  OVER  fiELflXflTIOM  FACTOR 


^  EQ.  4--7 
(  I  4  OJ  4  3-) 


( 


(n+  i)  tm 


flppfi  CXIO7flTI0N  OF 


D 


s 


*  denotes  the  n+1  th.  solution  using  Gauss-Seidel  without 
over-relaxation . 

The  iteration  process  has  converged  when  : 
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where  C  is 


S 


the  specified  error  for  all 


seam  elements . 


4.2  Multiple  Seams 

The  basic  method  is  the  same  as  for  a  single  seam,  except 
that  each  individual  seam  has  a  specified  depth,  thickness, 
material  properties  and  element  types  (see  Fig.  4.2) 
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For  M  seams  with  N  elements  in  each  seam,  the  stresses  at  any 
arbitrary  point  (x,y)  in  the  half-space  y^O  may  be  written  : 


%!*•*>  = 


•a 


G-  Cae-i-t-a  u")  -  (jCx-i 

V  J  14. 


n*i 


EQ.  4-.  8 


where  : 

J 

x  is  the  x  co-ord.  of  the  mid-point  of  the  j  th.  element  (all 
jfck  elements  have  the  same  x  co-ord.) 

m 

Functions  G  to  G  refer  to  functions  G  to  G  given  in  Section  3 

5  6  3  6 

and  f-ei.  1.  ,  for  the  m  th.  seam,  with  h  replaced  by  the  depth 

m 

of  the  seam  h  . 


The  solutions  are  obtained  by  solving  the  following  system 
of  equations  in  the  same  manner  as  for  a  single  seam  : 


£•<?.  I+.T 
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where 
seam 
i  .  e . 


I  ( 

K  and  K  are  the  shear  and  normal  stiffness  for  the  l  th. 


EQ.  10 


(A)  are  the  required  influence  coefficients  at  the  l  th.  seam, 
originating  from  the  m  th.  seam. 

(D)  and  (D)  are  the  shear  and  normal  displacement  discontinuities 
at  the  j  th.  element  of  the  m  th.  seam. 

This  is  a  system  of  2xNxM  linear  simultaneous  equations. 


The  program  DDSEAMS  has  been  developed  for  solution  of 
multiple  seam  mining  problems,  with  up  to  five  seams  and  80 
elements  modelling  each  seam. 
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Chapter  5 
PROGRAM  OPERATION 

The  computer  program  DDSEAMS  solves  the  numerical  problem  in 
Section  4.,  and  outputs  results  in  tabular  or  graphical  form.  The 
program  was  developed  to  be  simple  and  quick  to  set  up  a 
particular  mine  layout  and  allows  interactive  investigation  of 
the  stress  and  displacement  distributions  formed.  If  required  the 
program  may  be  run  in  batch  mode,  but  only  tabular  results  may  be 
printed  in  this  mode  as  the  graphics  option  must  be  run  from  a 
terminal.  A  list  of  the  major  subroutines  and  variables,  together 
with  a  program  listing  is  given  in  Appendix  1. 

The  program  is  written  in  Fortran  IV  and  is  intended  to  be 
used  in  conjunction  with  the  integrated  graphics  (*IG)  and 
plotting  (*PL0TLIB)  subroutine  packages  available  in  the 
University  of  Alberta  Computing  Centre.  Single  precision  is  used 
throughout,  as  this  gives  sufficient  accuracy  and  is  more 
efficient  in  computing  time  than  the  use  of  double  precision 
storage  and  calculations.  The  problem  solution  may  be  saved  as  a 
binary  data  file  for  re-running  the  program  again,  avoiding  the 
cost  of  re-solution. 

As  many  as  five  seams  with  up  to  eighty  elements  per  seam 
may  be  modelled  with  the  program  in  its  present  form.  Each  seam 
has  the  same  start  and  end  x  co-ords,  and  the  same  number  of 
elements.  The  equations  are  solved  using  Gauss-Seidel  iteration 
with  over-relaxation,  with  the  user  controlling  the  solution 
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tolerance,  over-relaxation  factor  and  the  maximum  number  of 
iterations.  After  the  solution  has  been  found,  and  the  output 
mode  has  been  defined,  there  are  four  alternative  procedures  : 
saving  the  solution  as  a  data  file,  calculating  the  stresses  and 
displacements  at  a  particular  seam,  caculating  the  stresses, 
strains  and  displacements  at  a  particular  depth,  and  the 
termination  of  the  program  run. 

5.1  Input/Output 
INPUT 

UNIT  4  ....For  a  new  run  the  data  defining  a  par t i cu 1 ar  prob lem 

is  read  in  from  unit  4. 

UN I  T  5  ....The  commands  controlling  the  program  run  are  read  in 

from  unit  5  which  defaults  to  ♦SOURCE*  (The  Terminal)  for 
interactive  operation. 

UNIT  7  ....For  a  re-run  the  solution  for  the  problem  is  read  in 

from  the  data  file  attached  to  unit  7. 

OUTPUT 

UNIT  3  ....The  offseam  results  are  written  to  the  file  attached 

to  unit  3  in  a  form  that  can  be  readily  used  in  any  later 
analyses  (e.g.  for  contour  plotting  of  principal  stress 
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calculations ) . 

UNIT  6  ....The  user  prompt  and  error  messages,  plus  graphical 

output  whilst  in  graphics  mode,  are  output  to  unit  6  which 
defaults  to  ♦SINK*  (The  Terminal)  for  interactive  use. 

U N If  7  ....  The  so  1 u  t i on  can  be  saved  on  unit  7  as  a  binary  data 

file. 

UNIT  8  ..The  printed  results  are  output  in  tabular  form  to  unit  8 
which  is  queued  for  printing,  a  typical  output  is  shown  in 
Appendix 

UNIT  9  ....The  data  for  plotting  of  graphical  output  is  written 

to  unit  9  whilst  in  interactive  graphics  mode. 

The  data  file  on  unit  4  must  contain  data  in  the  following 
form  : 

1  .  Title  Card  (20A4) 

A  title  consisting  of  less  than  80  alphanumeric  characters  (only 
the  first  35  characters  will  be  printed  in  graphical  mode). 

2.  Dimension  Card  ( 2 1 5 ,F 1 0 . 3 ) 

( i )  NSEAM  =  Number  of  seams 

(ii)  NSEG  =  Number  of  elements 

(iii)  HW  =  Half-width  of  each  seam  element 
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3.  Elastic  Constants  Card  ( 4E  1  1  . 4 , 3F5 . 2 ) 

(i)  EXX  =  Modulus  of  elasticity  in  x  direction 

(ii)  EYY  =  Modulus  of  elasticity  in  y  direction 

(iii)  EZZ  =  Modulus  of  elasticity  in  z  direction 

(iv)  GXY  =  Shear  Modulus  in  x,y  plane 

(v)  VXY  =  Poissons'  ratio,  y  strain  due  to  x  strain 

(vi )  VZX  =  Poissons'  ratio,  x  strain  due  to  z  strain 

(vii)  VYZ  =  Poissons'  ratio,  z  strain  due  to  y  strain 

4.  Rock  Parameters  Card  (2F10.3) 

(i)  DENS  =  Density  of  country  rock  (fo«.ce /uio it  volume) 

(ii)  SRATIO  =  Primitive  horizontal  to  vertical  stress  ratio. 

Seam  Cards  for  each  seam  : 


5.  Seam  Data  Card  ( F6 . 2 , F8 . 2 , 2E 1 1 . 4  ) 

(i)  THICK  =  Seam  thickness 

(ii)  DEPTH  =  Seam  depth 

(iii)  ESEAM  =  Stiffness  modulus  of  seam 

(iv)  GSEAM  =  Shear  modulus  of  seam 

6.  Seam  Excavation  Card  (8011) 

Each  character  represents  a  seam  element,  1=mined  and  2=unmined, 
NSEG  seam  elements  must  be  specified. 
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5.2  Running  Program 

In  interactive  mode  the  program  is  first  compiled  to  form 
the  object  code  for  the  program  : 


$r  *fortg  scards=ddseam. s  spunch=ddseam . o  par=d, noload 
The  program  can  now  be  run  as  follows  : 


$r  ddseam.o+*ig+*plot 1 ib  3=(results)  4=(data)  7=(save)  8=(print) 
9= (pdf ) 

The  program  prompts  the  user  for  execution  commands  and 
prints  the  alternatives  which  are  available  or  the  formats  of  the 
data  it  requires,  making  the  running  of  the  program  relatively 
s imp le . 


For  graphical  output,  after  graphical  data  has  been  written 
onto  unit  9  =  (pdf),  the  *calcompq  plot  program  must  be  run  in 
order  to  get  a  hard  copy  from  the  calcomp  plotter. 


Run  time  depends  upon  the  number  of  elements  and  seams,  and 
upon  the  number  of  internal  points  that  are  calculated.  For  four 
seams  with  75  elements  per  seam,  about  12  seconds  is  required, 
whereas  for  two  such  seams  only  6  seconds  is  required. 


References : 

1.  Crouch,  Steven  1.;  "Analysis  of  Stresses  and 

Displacements  Around  Underground  Excavations:  an 
Application  of  the  Displacement  Discontinuity  Method"; 
Dept •  of  Civil  and  Mineral  Engineering,  University  of 
M i nnesota ,  Minneapolis,  (1976). 
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APPENDIX  1. 

Computer  Program  DDSEAMS 

Internal  Subrout i nes 

ONSEAM. . .Computes  the  stresses  and  displacements  at  the  mid-point 
of  all  elements  in  a  given  range  for  a  given  seam. 

OFSEAM. . .Computes  the  stresses,  strains  and  displacements  for  a 
specified  depth  at  x  co-ords  the  same  as  the  centres  of  each  seam 


e  1  ement 

in  the  given  range. 

GFUNC . . 

.Computes  the  G  functions  as  given  in  ref.  1. 

GRAPH.  . 

.Plots  up  to  5  curves  on  one  graph 

GETTY.  . 

.Sets  up  the  arrays  for  subroutine  GRAPH 

STATS. . 

.Plots  the  problem  specifications 

MINEX.  . 

.Plots  a  cross  section  of  the  mine  layout 

I  BAR ...  Hatches  the  unmined  portions  of  the  seam  plotted  by  MINEX 


MESSGE. . .Places  a  message  on  the  plot 


. 
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Internal  Fund i ons 

VAVG. . .Calculates  an  average  value  for  a  real  array 
VMIN. .. Finds  a  minimum  vaue  for  a  real  array 
VMAX... Finds  a  maximum  value  for  a  real  array 

I TEM . . . P 1  aces  a  menu  item  at  (x,y)  and  returns  an  integer  picture 
name 

External  Subroutines 

PLOTS . . . Ini t i a  1 i ses  plotting  file 

1GCTRL ...  Device  control  routine 

1 GBGNS . . . Crea tes  a  sub-picture  within  an  existing  one  or  empty 
current  one 

IGENDS. . .Ends  a  sub-picture  and  returns  to  the  picture  started  by 
IGBGNS  or  IGCTNS 

IGVWPT ...  Viewport  transformation  for  a  picture  of  a  specific  size 
and  location 

IGCTNS ...  Re-act ivates  a  picture  to  the  current  one  for  additions 
or  changes 

IGVEC. . .Multiple  line  drawing  routine 

IGTXT ...  Draws  a  text  string  (transformation  of  picture  is 
appl i ed ) 


' 


168 


I GFMT . . . Conver t s  a  variable  vatue  to  a  character  representation  as 
text 

I GMR • • • Adds  a  non- i ntens i f i ed  1 i ne  to  the  end  of  the  current 

picture  (co-ord  displacement) 

IGTRAN ...  Linear  transformation  of  the  picture  when  drawn 

I GDRON . . . Tr ans 1  a tes  data  to  form  a  picture  at  the  terminal  or  to 
the  tape  unit  9  for  plot  files 

IGDELS. . .Deletes  a  picture  from  the  data  structure 

IGMA...Adds  a  non- intensi f ied  line  to  the  end  of  the  current 

picture/object 

IGPIKS. .. Returns  a  picture  name  for  the  current  position  of  the 
cursor  on  the  screen 

SCALE  ... F i nds  the  scale  parameters  for  graphs 

LINE... Plots  data  onto  graph  with  scaling  applied 

AXIS2 ...  Draws  a  single  graph  axis  with  marks  and  numeration 

IGTXTH . . . Adds  a  text  string  to  a  picture  without  transformations 
appl ied 

Main  Variables 

NSEAMS  -  total  number  of  seams 
NSEG  -  total  number  of  seam  elements 
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HW  -  half-width  of  each  seam  element 
DENSE  -  density  of  the  country  rock 

SRATIO  -  primitive  horizontal  to  vertical  stress  ratio 
DEPTH  -  array  containing  the  seam  depths 
THICK  -  array  containing  the  seam  thicknesses 

MIN  -  array  of  the  mined  state  for  each  seam  (1=mined,  2=unmined) 
ESEAM  -  elastic  moduli  for  each  seam 
GSEAM  -  shear  moduli  for  each  seam 
H  -  array  of  offseam  levels 

EXX.EYY.EZZ  -  elastic  moduli  for  country  rock  in  x,  y  and  z 
directions  respectively 

VYX.VZX.VYZ  -  Poissons'  ratios  for  the  country  rock 
DS  -  array  of  shear  displacement  discontinuities 
DN  -  array  of  normal  displacement  discontinuities 
ASS , ASN , ANS , ANN  -  Stress  Influence  Coefficients 
BSS , BSN , BNS , BNN  -  Displacement  Influence  Coefficients 
MENU1  -  array  of  1st.  command  menu  picture  names 
MENU2  -  array  of  2nd.  command  menu  picture  names 
IPLOT  -  picture  name  of  plotting  area  overlay 
ITEXT  -  picture  name  of  command  area  overlay 
GRAF  -  array  of  picture  names  of  plotted  graphs 
IMINE  -  picture  name  of  mine  cross  section 
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Program  DDSEAMS  listing 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c** 


c 


c 


TITLE 


'Analysis  of  Stresses  and  Displacements  Around  Multiple, 
Parallel  Seam-type  Deposits  by  the  Displacement 
Discontinuity  Method" 


CODE  NAME:  DDSEAM 
AUTHOR:  Steven  L.  Crouch 
IMPLEMENTED/ADAPTED  by:  Fred  E.  Eves 


VERSION:  O 


DATE:  FEBRUARY,  1978 


C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

COMMON/CONS/HW, DENSE , SRATIO , NSEAMS , NSEG 


Program  DDSEAM  computes  the  displacements  and  stresses  induced 
by  the  minig  of  multiple,  parallel  seam-type  deposits  using 
Crouch's  "Displacement  Discontinuity  Method".  The  program  can  be 
run  either  interactively  or  in  batch  mode.  The  computed  results 
may  be  obtained  in  the  form  of  printed  tables  and/or  plotted 
graphs.  The  *IG  and  *PLOTLIB  library  routines  are  used  to  imple¬ 
ment  the  graphics  option. 

Aside  from  the  changes  made  to  incorporate  the  plotting 
capability,  parts  of  the  rpogram  have  been  rewritten  inorder  to 
increase  efficiency. 


SUBROUTINES 
ONSEAM  - 
OFSEAM  - 

GFUNC 

GRAPH 

GETTY 

STATS 

MINEX 


called  by 
computes 
computes 
para  1  lei 
computes 


DDSEAM : 

the  stresses  &  displacements  along  a  seam, 
the  stresses  &  displacements  along  any 
line  in  the  half-space  y<0  not  on  a  seam, 
the  G-function  values  as  indicated  on  pages 
232-234  of  Crouch's  report. 

plots  up  to  5  curves  on  a  single  labeled  graph, 
moves  the  curve  data  into  columns  of  a  single  matrix 
for  the  curve  plotting  routine,  GRAPH, 
plots  the  specifications  describing  the  multiple 
seam  mine. 

plots  a  cross-section  of  the  multiple  seam  mine. 


COMMON/ CMO/ 
COMMON/CM  1/ 


1 


I  CARD (20) , EXX , EYY , EZZ , VYX , VZX , VYZ , GXY 
GAM  1 , GAM2 , 0 1 ,Q2,CN1 ,CN2,S1 1 , S22 , S 1 2 , C66 , B , BC66 , 

GM 1 2 , GM 1 2$2 , GM 1 $2 , GM2$2 , GM 1 $3 , GM2$3 , 0 1 P 1 ,Q2P1 
COMMON/ CM2/  ASS ( 80 , 5 , 5 ) , ASN( 80 , 5 , 5 ) , ANS(80, 5 , 5) , ANN (80, 5, 5) , 

BSS ( 80 ,5,5),BSN(80,5 , 5) , BNS ( 80 , 5 , 5 ) , BNN( 80 , 5 , 5 ) 
COMMON/ CM3/  DEPTH ( 5 ) . PS ( 5 ) , PN( 5 ) , THI CK ( 5 ) , E SE AM( 5 ) , 
GSEAM(5),MIN( 80 , 5  ) 

DS(80,5) ,DN(80,5) 

UXP(80,5) ,UXM( 80,5) ,UYP (80,5) ,UYM( 80,5) ,SIGS( 80,5) , 
SIGN( 80, 5 ) 

COMMON/ CM6/  UX ( 80 , 5 ) , UY ( 80 . 5 ) , ESXX(80,5) , E S Y Y ( 80 . 5 ) , E SXY ( 80 , 5 ) , 
SIGXX(80,5) ,SIGYY( 80 , 5 ) ,SIGXY(80,5) ,H(5) 


COMMON/ CM4/ 
COMMON/ CM5/ 


DIMENSION  X(80) ,XX(82) , YY(248) 

COMPLEX  *8  Y Y I D ( 3 ) , CDUMMY 

REAL  K 1 ,K2,KS,KN,XVWPT(5), Y VWPT ( 5 ) 

INTEGER  WINDOW, CALC .DAT SET , I DAT ( 4 ) , MENU  1 ( 5 ) ,MENU2(5)  , 


' 
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1  GRAF(3) , XTITLE(5) , YTITLE 

DATA  I  DAT/ 'SAVE' ,  'ONSE'  ,  'OFFS'  ,  'STOP'/, 

1  ITEXT/'TEXT'/,  I  PLOT/ ' DRAW ' / ,  I DENT/ ' I DEN ' / , 

2  IMINE/'MINE'/,  I  ST AT/ ' STAT ' / ,  CALC/ ' CALC '/ , 

3  GRAF/'PLT 1 ' , 'PLT2' , ' PLT3 ' / , 

4  WINDOW/'WIND'/, 

5  XVWPT/-1 .  ,-1  .  ,  .5,  .5,-1 ./,  YVWPT/- .75,  .75,  .75, -.75, -.75/, 

6  XT  I TLE/ ' HOR I '  ,  ' ZONT '  ,  ' A  L  D '  ,  ' I  ST  A '  ,  ' NCE ' / ,  YTITLE/'  '/ 


-  initialize 

LVL  =0 


Determine  whether  a  new  problem  is  to  be  solved  or  whether  a 
previously  solved  problem  is  to  be  read  in. 

I0PT=1  ->  pre-defined  problem. 

I0PT=2  ->  new  problem. 

5  PRINT  2001 

READ  101 1 ,  I  OPT 
GOTO  (10, 20 ) , I OPT 
PR  I  NT  2096 ,  I  OPT 
GOTO  5 

10  READ( 7 )  ASS, ANN, ASN.ANS, BSS.BNN.BSN.BNS, DS,DN, PS, PN,X,H, DEPTH, 

1  THICK, ESEAM.GSEAM, EXX , EYY , EZZ, GXY ,VYX,VZX,VYZ,K1 ,K2, 

2  GAM  1 , GAM2 , Q 1 ,Q2,CN1 .CN2.S1 1 , S22 , S 1 2 , C66 , B , BC66 , 

3  GM 1 2 , GM 1 2$2 , GM 1 $2 , GM2$2 , GM 1 $3 , GM2$3 , 0 1 P 1 ,Q2P1 , 

4  DENSE, SRATIO.HW.DMAX, OMEGA, 

5  MIN, ICARD .NSEAMS ,NSEG, NITER 

allow  for  a  change  in  the  pre-existing  stress  ratio. 


PRINT  2003 

READ  1002 ,  SRATIO 

GOTO  500 


-  Read  in  parameters  defining  a  new  problem. 

20  READ( 4 , 1001 )  ICARD 

READ( 4 , 1003 )  NSEAMS , NSEG , HW 

READ( 4 , 1004)  EXX , EYY , EZZ , GXY , VYX , VZX , VYZ 

READ( 4 , 1005 )  DENSE, SRATIO 

-  Compute  constants  GAM1,  GAM2,  etc. 

P I =4 . *ATAN( 1  .  ) 

SI  1  =  1  ./EXX 
S22  = 1 . /EYY 
S33= 1 . /EZZ 
S 1 2=-VYX/EXX 
S  1  3  =  -VZX/EXX 
S23=-VYZ/EZZ 

C0N=S1 1*S22-S12*S12+(2 . *S 1 2*S 1 3*S23-S 1 1 *S23*S23-S22*S 1 3*S 1 3 )/S33 
Cl 1=(S22-S23*S23/S33)/C0N 
C 12=-(S 12-S13*S23/S33 )/CON 
C22=( SI 1-S13*S13/S33)/C0N 
C66=GX Y 

DET=C 1 1*C22-C12*C12 
IF  (DET .GT .0. )  GO  TO  50 
PR  I  NT 20 1 7 
STOP 

50  K  1  =  SORT ( C22/C 1 1 ) 

K2  =  0 . 5* ( C 1 1*C22-C12*(2. *C66  +  C12) )/ (Cl 1*C66) 

IF  (K2-K1)  55,60,65 
55  PR I NT 20 1 8 


' 
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STOP 

60  PRINT 20 1 9 
STOP 

Compute  constants  used  in  subroutine  GFUNC 
65  SUM=SQRT( (K1+K2)/2.  ) 

DI F  =  SORT ( ( K2-K 1  )/2 .  ) 

GAM  1 =SUM+DI F 

GAM2=SUM-DIF 

GM 12=GAM1*GAM2 

GM 1 2$2=GM 12**2 

GM 1 $2  =  GAM 1**2 

GM2$2=GAM2**2 

GM 1 $3  =  GM 1 $2  *GAM 1 

GM2$3  =  GM2$2  *GAM2 

01= (Cl 1*GM1$2-C66)/ (C12+C66) 

02= ( C 1 1 *GM2$2-C66 )/ (C12+C66) 

01P 1 = 1 . +01 
02P 1 = 1 . +Q2 

B  =  Q 1 P 1 *Q2P  1  / ( 2 . *PI*(Q1-Q2) ) 

CN1 =(GAM1+GAM2)/(GAM1 -GAM2 ) 

CN2  =  2 . *GM 12/ (GAM1-GAM2) 

BC66  =  B  *C66 


-  Read  in  specifications  for  the  various  seams. 

DO  70  M= 1 , NSEAMS 

READ ( 4 , 1006 )  THI CK ( M ) , DE PTH ( M ) , ESEAM(M) ,GSEAM(M) 

READ( 4, 1011)  (MIN(N.M) ,N= 1 ,NSEG) 

70  CONTINUE 

******************************************************* 

***  COMPUTE  INFLUENCE  COEFFICIENTS.  *** 

******************************************************* 

DO  200  L  = 1 .NSEAMS 
HL=DEPTH( L ) 

DO  200  M= 1 .NSEAMS 
HM=DEPTH( M ) 

DO  200  K= 1 , NSEG 
KK  =  2* (K- 1  ) 

XP=HW* ( KK+ 1 ) 

XM=HW* ( KK- 1 ) 

CALL  GFUNC ( XP . HM , -HL , G 1 P . G2P , G3P , G4P , G5P , G6P , G7P , G8P , G9P , 

*  G10P) 
CALL  GFUNC (XM.HM, -HL , G 1 M . G2M , G3M , G4M , G5M , G6M , G7M , G8M , G9M , 

*  G 1 OM ) 
ASS(K,M,L)=BC66*(G5P-G5M) 

ASN(K, M, L )=BC66*(G6P-G6M) 

ANS(K,M, L) =BC66  * (G3P-G3M) 

ANN( K,M,L)=BC66*(G4P-G4M) 

BSS(K,M. L)=B*(G7P-G7M) 

BSN(K,M, L)=B*(G8P-G8M) 

BNS(K,M, L)=B*( G9P-G9M ) 

BNN(K,M, L)=B*(G10P-G10M) 

200  CONTINUE 

INITIALIZE 

210  DO  220  M= 1 .NSEAMS 
PS ( M ) =0 . 

PN(M)=DENSE*DEPTH(M) 
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DO  220  1=1 , NSEG 
DS ( I , M ) =0 . 

DN( I , M ) =0 . 


220  CONTINUE 
C 


230  DO  240  1=1 , NSEG 


X( I )=HW*(2*I-1-NSEG) 


240  CONTINUE 

C*************************** 


************************** 


***************** 


C*** 

C*** 


SOLVE  EQUATIONS  BY  GAUSS-SIEDEL  ITERATION 
PROCESS  -  WITH  OVER-RELAXATION. 


********************** 


solve  for  shear  &  normal  displacement  d i cont i nu i t i es . 

NITER=0 

PRINT  2035 

250  PR  I  NT 2042 

READ  1013,  N.TOL, OMEGA 

251  PRINT  2043 
READ  1000,  IANS 
GOTO( 252 , 254 ) , IANS 

GOTO  251 

252  PR  I NT 2044 

254  DO  290  ITER=1,N 
DMAX  =0 . 

DO  280  L= 1 .NSEAMS 


HO=THI CK( L ) 

KS=GSEAM(L)/HO 
KN=ESEAM(L)/HO 
PPS  =  PS( L  ) 

PPN=PN( L ) 

ASSLL  =  ASS(  1 ,L,L) 

ANNLL=ANN( 1 , L , L ) 

DO  280  I L= 1 , NSEG 
MINE  =  MIN(IL, L  ) 

TMS=0 . 

TMN=0 . 

DDS=DS ( I L , L ) 

DDN=DN ( I L , L ) 

DO  260  M= 1 .NSEAMS 
DO  260  JM= 1 .NSEG 
DMS=DS(UM,M) 

DMN=DN( JM , M ) 

LD= I L-JM 
K= I ABS ( LD )  +  1 
F  AC= I S I GN( 1 , LD) 

TMS=TMS  -  ASS ( K . M , L ) *DMS  -  FAC*ASN(K,M,L) *DMN 
TMN=TMN  -  FAC*ANS(K,M, L ) *DMS  -  ANN( K , M , L ) *DMN 


260 


CONTINUE 


C 


265 


GO  TO  (265,270) .MINE 

TMS=TMS-PPS 

TMN=TMN-PPN 

DELS=TMS/ASSLL 

DELN=TMN/ANNLL 


C 


check  8>  correct  for  complete  closure  boundary  condition. 
IF  (DDN.LT.HO.OR.TMN.GT.O. )  GOTO  275 
DDN=HO 
DELN=0 . 

GOTO  275 


C 


■ 


on  non 
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270  TMS=TMS  +  KS*DDS 

TMN=TMN  +  KN*DDN 
DELS=TMS/( ASSLL  -  KS) 
DELN=TMN/ (ANNLL  -  KN) 


275 


C 


280 

282 

284 

290 


DS ( I L , L  )  =DDS  +  OMEGA  +  DELS 
DN(IL,L)=DDN  +  OMEGA  *DE  LN 

DELS=ABS(DELS) 

DELN=ABS( DELN ) 

DMAX= AMAX 1 ( DMAX , DELS, DELN) 
CONTINUE 

GOTO  ( 282 , 284 ), IANS 
PR I NT  2045 ,  ITER, DMAX 
IF  (DMAX.LE.TOL)  GOTO  295 
CONTINUE 


NI TER=NI TER+N 
PRINT  2039 ,  NITER, DMAX 
292  PR  I  NT 2040 

READ  1000 ,  IANS 
GOTO ( 250, 500, 9999) , IANS 
PRINT  2096 ,  IANS 
GOTO  292 


295  NITER=NITER+ITER 

PR I  NT  204 1 ,  NITER, OMEGA, DMAX 
C 

C**********************************************^t***t***********i:****** 

c***  *** 

c***  COMPUTATION  CONTROL  SECTION  *** 

C***  *** 

c********************************************************************** 
C  Determine  mode  of  output. 

C  M0DE=1  ->  graphical  only. 

C  M0DE  =  2  ->  both  graphical  8>  printed. 

C  M0DE=3  ->  printed  only. 

500  PR I NT 2070 

READ  1000, MODE 
GOTO  (505,505,505),  MODE 
PR  I NT 2096 ,  MODE 
GOTO  500 

505  PR  I  NT 207 1 ,  MODE 

GOTO  (515,510,510),  MODE 


OUTPUT  THE  MINE  SPECS. 


510  WR I TE ( 8 , 2002 )  ICARD 

WR I TE ( 8 , 2008 )  NSE AMS , NSEG , HW 
WR I TE ( 8 , 2009 )  DENSE, SRATIO 
WR I TE ( 8 , 2010 )  EXX,EYY,EZZ,GXY, VYX , VZX 
WRITE (8, 2021 )  K1 ,K2,GAM1 ,GAM2,Q1 ,02 
DO  511  M= 1 .NSEAMS 

WRITE (8 , 2012 )  M, THICK) M) , DEPTH( M ) 
WRITE (8, 2013)  (MIN(N,M)  , N= 1 , NSEG ) 

511  CONTINUE 

WRITE(8 , 204 1 )  NI TER , OMEGA , DMAX 
IF  ( MODE . EO . 3 )  GOTO  525 


VYZ 


ESEAM(M) , GSEAM(M) 


***  SET  UP  GRAPHICAL  PICTURES 


o  o  o 
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c 


515  CALL  PLOTS 

CALL  IGCTRL(CALC, 'PACK' , 
CALL  IGCTRL( 'TERM' , 'KEEP 
YTOP= . 75 

CALL  IGBGNS( ITEXT ) 

C -  menu  no.  1 

MENU  1 ( 1 )  = 

MENU  1 ( 2 )  = 

MENU  1(3)=  ITEM(  ' 

MENU  1(4)=  ITEM(  ' 

C -  menu  no.  2 

MENU2 ( 1 ) =  ITEM( ' 

MENU2 ( 2 ) =  I TEM ( ' 

CALL  IGENDS( ITEXT) 

CALL  IGVWPT( ITEXT , 

CALL  I GBGNS ( I  PLOT  ) 

C -  define  mine  cross-section 


'ON'  ) 
'  .  1  ) 


ITEM( 'BLOW-UP  GRAPH<E> ' , 
ITEM( 'CALCOMP  COPY<E>',. 
'MINE  X-SECTION<E> ' 


.55, YTOP- .03) 
55, YTOP- .09) 

,  .55, YTOP- . 15) 


CONT I NUE  <E  > ' ,  . 55, YTOP- .21) 

REDRAW<E> ' , .05, YTOP- .03) 
CALCOMP  COPY<E>' , .05, YTOP- .09) 

.5,  1  .  ,-.75,  .75) 


CALL  MINEX(IMINE) 

CALL  IGVWPT( IMINE , -  1 .  ,  1 .  , 1 .  , 3 .  ) 

C -  define  statistics  display. 

CALL  ST  ATS ( I  ST  AT  ) 

CALL  I GCTNS ( I  ST  AT ) 

CALL  IGBGNS( IDENT ) 

CALL  IGENDS( IDENT) 

CALL  I GENDS ( I  ST  AT ) 

CALL  I GVWPT ( ISTAT.O.O, 1.0, 0.0, 1.0) 

C -  define  viewports  for  each  of  the  graphs. 

CALL  I GBGNS ( GRAF ( 1  )  ) 

CALL  I GENDS ( GR A F (  1  )  ) 

CALL  I GVWPT ( GRAF ( 1), "I, 0.0, 0.0, 1.0) 
CALL  I GBGNS (GRAF(2) ) 

CALL  I GENDS (GRAF(2)  ) 

CALL  I GVWPT (GRAF (2 ) , -  1  .  ,0. 0,-1. 0,0.0) 
CALL  I GBGNS (GRAF(3)) 

CALL  I GENDS ( GRAF(3)  ) 

CALL  I GVWPT (GRAF(3) ,0.0, 1.0, -1.0, 0.0) 
CALL  I GENDS ( I  PLOT  ) 


-  place  drawing  area  in  largest  left  square  of  screen. 

-  also  box  this  region. 

CALL  I GVWPT ( I  PLOT , -  1  . 0 , O . 5 , - . 75 ,  . 75 ) 

CALL  I GVEC ( 5 , XVWPT , YVWPT  ) 

525  PR  I NT 2030 

READ  1 008 ,  DATSET 
DO  550  1=1,4 

IF  (DATSET . EO. IDAT( I  )  )  GOTO  575 
550  CONTINUE 

PRINT  2000,  DATSET 
GOTO  525 

575  GOTO  (600,700,800,9999) , I 
C*** 

C***  SAVE  the  computations  of  influence  coefficients,  displacement 
C***  discontinuity,  and  other  variables  for  another  run. 

C  *  *  * 

600  WR I TE ( 7  )  ASS, ANN, ASN , ANS , BSS , BNN , BSN.BNS, DS.DN.PS, PN,X,H, DEPTH, 

1  THICK,ESEAM,GSEAM,EXX,EYY,EZZ,GXY,VYX,VZX,VYZ,K1 , K2 , 

2  GAM  1 , GAM2 , 0 1 .Q2.CN1 , CN2 , S 1 1 , S22 , S 1 2 , C66 , B , BC66 , 

3  GM 1 2 , GM 1 2$2 , GM 1 $2 , GM2$2 , GM 1 $3 , GM2$3 , 0 1 P 1 , Q2P1 , 

4  DENSE, SRATIO.HW.DMAX, OMEGA, 


■ 


oo  on  on  non 
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5  MIN, ICARD.NSEAMS.NSEG, NITER 

GOTO  525 

C  *  *  * 

C***  COMPUTE  stresses  and  displacements  at  seam  level (s). 

C*  *  * 

700  PRINT  2080,  NSE AMS , N5EG 
READ  1009 ,  ISM, IB, IE 

IF  ( ISM . GE . 1 . AND . ISM. LE .NSEAMS)  GOTO  710 
PR I  NT  208 1 ,  ISM 
GOTO  700 

710  IF  (IB.GE.1  .AND.  I E . GE . 1  .AND.  IB.LE.NSEG  .AND.  IE.LE.NSEG 
1  . AND .  IB . LE . IE )  GOTO  720 

PRINT  2095 ,  IB, IE 
GOTO  700 
C 

720  CALL  ONSEAM( ISM, IB , IE ) 

GOTO  (775,725,725),  MODE 
725  WRITE(8 . 2050)  ISM 
DO  750  U= I B , IE 

WRITE (8 , 2051 )  U,UXP(d, ISM) ,UXM( J, I SM ) , DS ( U , I SM ) , UYP ( J , I SM ) , 

1  UYM( J, ISM) ,DN(U, ISM) , SIGS( J, ISM) , SIGN( J , ISM) 

750  CONTINUE 
C 

775  GOTO  (776,776,525),  MODE 
GRAPHICAL  OUTPUT 

776  CALL  IGCTRL( 'TERM' , 'ERASE' ) 

CALL  I GBGNS ( I DENT ) 

CALL  I GMR ( .015, .02) 

CALL  I GTXT ( '<RSCL>' ,1.5, 'ONSEAM:<E>' ) 

CALL  IGTXT( '<RSCL>' , .6667, '  SEAM  <E>') 

CALL  IGFMT( ISM, ' I ' , 1 ) 

CALL  I GTXT ( '  SEGS  <E>' ) 

CALL  IGFMT ( IB , ' I ' , 2 ) 

CALL  I GTXT ( ' -<E> ' ) 

CALL  I GFMT ( I E , ' I ' ,2) 

CALL  IGENDS( IDENT ) 

N  =  0 

DO  780  <J=  I B  ,  IE 

N  =  N+  1 

XX(N)=X( J) 

780  CONTINUE 
N2=N+2 

plot  shear  displacements  and  displacement  discontinuity  values. 
CALL  GET Y Y ( 3 , UXP ( I B , ISM) , ' +U<BSUB>S<E> ' ,UXM( IB, ISM) , 

1  ' -U<BSUB>S<E>'  ,DS( IB , ISM)  ,  '  D<BSUB>S<E> ' , N2 , YY , Y Y I D ) 

CALL  GRAPH ( XX,YY,N2,3,YYID,XTITLE,  19.YTITLE,  1  , 

1  'SHEAR  DISPL.  &  DISPL.  DI SCONTINUITY<E> '  , GRAF ( 1  )  ) 

plot  normal  displacements  and  displacement  discontinuity  values. 
CALL  GETY Y ( 3 , UYP ( IB, ISM) , ' +U<BSUB>N<E> ' ,UYM( IB, ISM) , 

1  ' -U<BSUB>N<E>'  ,DN( IB , ISM) ,  '  D<BSUB>N<E> '  , N2 , YY , YYID  ) 

CALL  GRAPH(XX,YY,N2,3,YYID,XTITLE,  19.YTITLE,  1  , 

1  'NORMAL  DISPL.  &  DISPL.  D I SCONT I NU I T Y<E> '  , GRAF ( 2  )  ) 

plot  shear  and  normal  stresses  on  each  seam  element. 

CALL  GETYY(2,SIGS(IB,ISM), ' SIG<BSUB>XY<E> ' , SIGN( IB , ISM) , 

1  'SIG<BSUB>YY<E>' , DUMMY , CDUMMY , N2 , Y Y , Y Y I D ) 


■ 


on  on  o  o  o 
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CALL  GRAPH(XX,YY,N2,2,YYID,XTITLE, 1 9 , YTI TLE  , 1 , 

1  'SHEAR  &  NORMAL  STRESSES<E> ' . GRAF ( 3 ) ) 

GOTO  9000 

C*** 

C***  COMPUTE  stresses  and  displacements  at  specified  off-seam  locations 

C*** 

800  PRINT  2090,  NSEG 
LVL  =  L VL  + 1 

READ  1010,  H( LVL), IB, IE 
IF  (H(LVL).LE.O.O)  GOTO  810 
PRINT  209 1 ,  H(LVL) 

GOTO  800 

810  IF  (IB.GE.1  .AND.  I E . GE . 1  .AND.  IB. LE. NSEG  .AND.  IE. LE. NSEG 
1  . AND .  IB . LE . IE )  GOTO  820 

PRINT  2095 ,  IB, IE 
GOTO  800 
C 

820  CALL  OFSEAM( LVL , IB , IE ) 

GOTO  (875,825,825),  MODE 
825  WR I T E ( 8 , 2060 )  H(LVL) 

DO  850  J= IB , IE 

WRITE (8, 2061 )  U,UX(U, LVL ) ,UY(U, LVL) ,ESXX( J, LVL) , ESYY( J, LVL) , 

1  ESXY( J, LVL ) , SIGXX( J, LVL ) , SIGYY ( J , LVL ) , SIGXY ( J , LVL ) 

WRITE (3, 3061 )  U,H(LVL  )  , X ( J  )  , UX ( J , LVL ) ,UY( J, LVL  )  , ESXX( J, LVL) , 

1  ESYY(U, LVL) , ESXY( J, LVL) , S I GXX ( J , L VL ) , S I GY Y ( U , LVL ) , S I GX Y ( U , LVL ) 

3061  FORMAT  ( I  5 , 2F 10 .  1 , 5F 10 . 6 , 3F 10 .  1 ) 

850  CONTINUE 
C 

875  GOTO  (876,876,525),  MODE 
GRAPHICAL  OUTPUT 

876  CALL  IGCTRL( ' TERM' ,' ERASE ' ) 

CALL  I GBGNS ( I DENT ) 

CALL  I GMR ( .015, .02) 

CALL  I GTXT ( '<RSCL>' , 1 . 5 , ' OFFSEAM : <E> ' ) 

CALL  IGTXT(  '  <RSCL> '  ,  . 6667 ,  '  Y  =  <E>') 

CALL  IGFMT(H(LVL), ' F ' ,7, 1) 

CALL  I GTXT ( '  SEGS  <E>' ) 

CALL  I GFMT (IB, ' I ' , 2 ) 

CALL  I GTXT ( ' -<E > '  ) 

CALL  I GFMT ( I E , ' I ' , 2) 

CALL  I GENDS ( I DENT ) 

N  =  0 

DO  880  J= IB , IE 
N  =  N+  1 

XX  ( N )  =X  ( <J  ) 

880  CONTINUE 
N2=N+2 

plot  strain  at  this  level. 

CALL  GETYY(3, ESXX( IB , LVL) , ' E<BSUB>XX<E> ' , ESYY ( I B , LVL ) , 

1  ' E<BSUB>YY<E> '  , E SX Y ( I B , LVL  )  ,  ' E<BSUB>XY<E> '  ,N2, YY.YYID) 

CALL  GRAPH( XX,YY,N2,3,YYID,XTITLE, 19 , YTITLE , 1 , 

1  'STRAIN  S<E> ' , GRAF (  1  )  ) 

plot  stresses  at  this  level. 

CALL  GETYY(3, SIGXX( IB, LVL) , ' SIG<BSUB>XX<E> ' , S I GYY ( I B , LVL ) , 

1  'SIG<BSUB>YY<E>' ,SIGXY(IB,LVL), 'SIG<BSUB>XY<E>' ,N2, YY.YYID) 

CALL  GRAPH ( XX,YY,N2,3,YYID,XTITLE,  19, YTITLE,  1, 


oo  oo  o  onono  no 
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'S  T  R  E  S  S  E  S<E> ' , GRAF ( 2  )  ) 


P 


1  ot 

1 

1 


shear 

CALL 

CALL 


and  normal  displacements  at  this  level. 

GETVY ( 2 , UX ( IB , LVL  )  ,  ' U<BSUB>X<E> ' , UY( IB , LVL ) , 

' U<BSUB>Y<E>  '  , DUMMY , CDUMMY , N2 , Y Y , Y  Y I D ) 
GRAPH(XX,YY,N2,2,YYID,XTITLE,  19 , YTITLE ,  1, 

'SHEAR  &  NORMAL  D I SPLACEMENTS<E> ' , GRAF 


GOTO  9000 


(3) 


) 


***  GRAPHICS  MONITOR  *** 


8000  CALL  IGTRAN( ITEXT , 'WIND' , 0 . 0 , 0 . 5 , - . 75 , .75) 
CALL  IGCTRL( 'TERM' ,' ERAS' ) 

8010  I NDEX  =  IGPIKS(2, MENU 2 ) 

GOTO  (8030,8020) , INDEX 

8020  CALL  IGTRAN( '*MP*' , 'WIND' ,-1 . , .5, -.75, .75) 
CALL  I GDRON( 'CALC ' ) 

CALL  IGTRAN( ' *MP* ', 'SCALE 1 . ) 

GOTO  8010 

8030  CALL  IGTRAN( IPLOT , 'SCALE' , 1 .0) 

CALL  IGCTRL( 'TERM'  ,  'ERAS'  ) 

9000  CALL  1GTRAN( ITEXT , WINDOW ,0.5, 1 . O , - . 75 , .75) 

9010  I NDEX  =  IGPIKS(4 , MENU  1 ) 

G0T0(91OO, 9200, 9300, 9400) , INDEX 


9100  CALL  MESSGE(  '>>INDICATE  PLOT?<E> '  ,  .  5 1 , - YTOP ) 
I NDEX  =  IGPIKS(3,GRAF ) 

CALL  IGDELS( 'MESS' ) 

GOTO  (91 10,9120,9130) , INDEX 
9110  CALL  IGTRAN( IPLOT , WINDOW, -1 .0,0. 0,0.0, 1 .0) 
GOTO  8000 

9120  CALL  IGTRAN( IPLOT .WINDOW, -  1 .0,0.0, -  1  .0,0.0) 
GOTO  8000 

9130  CALL  IGTRAN( IPLOT .WINDOW, 0.0,  1  .0, -  1  .0,0.0) 


9200  INT=IGBGNS(0) 

CALL  I GMA (  .  5  1  , - . 67 ) 

CALL  IGTXT( ' SELECT  OPTION:<E>') 

NP I C 1  =  I T  EM( 'COMPLETE<CRLF>DISPLAY<E> ' ,  . 51 , - . 71 ) 
NP IC2=ITEM(  ' GRAPHS<CRLF>ONLY<E> '  ,  . 81  , - . 71 ) 

I NDEX  =  IGPIKS(NPIC1  , NP I C2 ) 


CALL 

I GENDS ( I NT ) 

CALL 

I GDE  LS ( I NT ) 

CALL 

I GTRAN(  ' *MP  * '  .WINDOW, - 

1 . , .5, -.75, 

GOTO 

(9210,9220),  INDEX 

9210  CALL 

I GDRON( CALC ) 

GOTO 

9230 

9220  CALL 

I GTRAN( IPLOT, WINDOW, -  1 

. ,0. ,0. . 1 . ) 

CALL 

I GDRON (CALC) 

CALL 

IGTRAN( I  PLOT, WINDOW, -1 

. ,0. , -1 . ,0. 

CALL 

IGDRON( CALC ) 

CALL 

I GTR AN ( I  PLOT , W I NDOW , 0 . 

.  1  .  , -  1  .  ,0.  ) 

CALL 

IGDRON(CALC) 

■ 


* 
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c 


CALL 


9230  CALL 
CALL 
GOTO 
C 
C 


9300  CALL 
GOTO 
C 
C 


9400  CALL 
GOTO 


C 


IGTRAN( I  PLOT,  'SCAL'  ,  1  .  ) 

I GTR AN (  ' *MP  *  '  ,  ' SCAL '  ,  1  .  ) 
I GCTRL (  'TERM'  ,  ' ERASE '  ) 
9010 


IGTRAN( I  PLOT , WINDOW, -  1  .  ,  1 .  ,  1 .  , 3 .  ) 
8000 


I GCTRL ( 'TERM' , 'ERASE' ) 
525 


9999  IF  ( MODE . LE . 2  )  CALL  I GCTRL ( CALC ,' ENDPLOT ' ) 
STOP 


C 

C 

C 

1000 

1001 

1002 

1003 

1004 

1005 

1006 

1007 

1008 

1009 

1010 
101  1 
1013 
C 

2000 

2001 

2002 

2003 

2008 

2009 

2010 


2012 


2013 

2017 

2018 
2019 
202  1 


FORMAT  STATEMENTS. 


FORMAT ( I  1 ) 

FORMAT  ( 20A4 ) 

FORMAT  (F10.0) 

FORMAT  (215, F 10.3) 

FORMAT  ( 4  E 1  1 .4.3F5.3) 
FORMAT  (2F10.3) 

FORMAT  (F6.2,F8.2,2E11 .4) 
FORMAT  (F8.2.2I4) 

FORMAT ( A4 ) 

FORMAT (315) 

FORMAT ( F 10.  1 ,215) 

FORMAT  (8011) 

FORMAT( 13, 2F 10.0) 


FORMAT ( ' 
FORMAT ( ' 
FORMAT ( 1 
FORMAT  ( 
FORMAT  ( 
1 ' IN  EACH 
FORMAT ( 

1 

FORMAT 

1 

2 

FORMAT 

1 

2 

3 

4 

FORMAT 
FORMAT 
1  .  '  ) 
FORMAT 
FORMAT 
FORMAT 

1 

2 

3 

4 

5 


DATASET  NAME  -  ' 
PREDEFINED?  (Y=1 


,A4//) 
N  =  2  )  . 


(II)') 


-***  BAD 
PROBLEM 
HI , 20A4 ) 

'OInput  new  SRATIO,  if  different.  ...(F10)') 

'ONUMBER  OF  SEPARATE  SEAMS  = '  , I  3 , / ' ONUMBER  OF  SEGMENTS ‘ 
SEAM  =' , I4,/'0HALF-WIDTH  OF  SEGMENTS  =',F6.2) 

OUNIT  WEIGHT  OF  ROCK  =',F6.3/, 

ORATIO  OF  PRE-EXISTING  STRESSES  SIGXX/SIGYY  =' 


'OELASTIC  CONSTANTS  EXX  = ' . E 1 1 
Ell. 4,/'  EZZ  = ' , E 1 1 . 4 , / '  GXY  =',E11 
F6 . 2 , / '  VZX  =  ' , F6 . 2 , / '  VYZ  =',F6.2) 
'-SPECIFICATIONS  FOR  SEAM  NUMBER', 12,/ 


. F6 . 2  ) 
EYY  = 
VYX  = 


'0 


(  'O 


SEAM  THICKNESS  =',F9.2/ 

DEPTH  = ' , F9 . 2/ 

MODULUS  OF  ELASTICITY  =',E10.3/ 

SHEAR  MODULUS  =',E10.3) 

MINING  PATTERN  (1  =  MINED,  2  =  UNMINED) 


,//5X,  1331 1  ) 

INADMISSIBLE  ELASTIC  CONSTANTS,  C  1  1  *C22 -C  1  2 *C  1  2  LE 


'-***  COMPLEX  GAM  1 ,  GAM2  NOT  ALLOWED 
'-***  GAM  1  =  GAM2  NOT  ALLOWED.') 
'OELASTIC  PARAMETERS  -  K1  =',F6.2,/ 
22X 
20X 
20X 
22X 
22X 


' K2  = ' , F6 . 2/ 


'  GAM  1 
'  GAM2 
'01  = 
'02  = 


= ' , F6 . 2/ , 
= ' , F6 . 2/ , 
, F6 . 2/ , 

, F6 . 2  ) 


O 


. 


180 


2030 

2035 

2039 


2040 


204  1 


2051 

2060 


2061 

2070 


207  1 
2080 

208  1 

2090 

2091 

2095 

2096 
C 

c 


FORMAT (  'OCOMMAND :  SAVE , ONSEam , OF FSeam , STOP .  .  . ( A4  )  '  ) 

FORMAT ( ' -  SOLVE  for  displacement  discontinuities.'/) 

FORMAT  ('0***  ITERATION  PROCESS  HAS  FAILED  TO  CONVERGE  ADEQUATELY 
1/'  AFTER  ',13,'  ITERATIONS.'/'  MAX  I MUM  D I F F E RENCE ' , 

2'  BETWEEN  SUCCESSIVE'/'  ITERATES  IS',E10.3) 

FORMAT( 'OINDICATE  ACTION:'/ 

1'  1  -  Continue  iteration'/ 

2'  2  -  Accept  solution  &  proceed'/ 

3'  3  -  Stop'/) 

FORMAT( ' -SOLUTION  DATA  FOR  GAUSS-SIEDEL  ITERATION  PROCESS'/, 

*  'O  NUMBER  OF  ITERATIONS :', I  10/ , 

*  '  RELAXATION  F ACTOR : ' , 3X , F 1 0 . 2/ , 

*  '  ATTAINED  TOLERANCE: ' ,2X,E10. 3/) 

#  i terat ions , tolerance .relax,  factor:  .  .  . ( 1 3 , 2F 10 . 0) ' / 

*  '&?') 

FORMAT( 'OPRINT  INTERMEDIATE  ITERATION  VALUES?  (1=Y,2=N)') 

FORMAT ( ' -  ITERATION  MAX.  RELAXATION'/) 

FORMAT ( 4X , I4.8X.E10.3) 

FORMAT  ('1  STRESSES,  DISPLACEMENTS 
1 TY  COMPONENTS  AT  SEAM  #',I1,/ 

2 ' 0  SEG  UX(POS)  UX(NEG)  DX  (=DS) 

3  SIGXY  SIGYY',/) 

FORMAT  (I5,6F10.6,2F10. 1 ) 

STRESSES  AND  DISPLACEMENTS  AT  OFF-SEAM  DEPTH:  Y=', 
SEG  X-DISP  Y-DISP  XX-STRAIN  YY-STRAIN  XY-STRAIN 
SIGYY  SIGXY'/) 


2042  FORMAT ( 


2043 

2044 

2045 
2050 


AND  DISPLACEMENT  DISCONTINUI 


UY(POS)  UY(NEG)  DY  ( =DN ) 


FORMAT  ( ' 1 

1  F8  .  1  // ' 

2  SIGXX 
FORMAT 
FORMAT 

1 
2 

3 

FORMAT 
FORMAT 
FORMAT ( 
FORMAT 
FORMAT ( 
FORMAT ( 
FORMAT ( 


END 


( I  5 , 5F 10 . 6 , 3F  10  .  1  ) 

('  INDICATE  OUTPUT  MODE  .  .  .  (  I  1  )  '  / 
'  1  =  graphic  only'/ 

'  2  =  graphic  &  print'/ 

'  3  =  print  only' ) 

( '  OUTPUT  MODE  '  ,  I  1  ) 


( 


INPUT 
*  *  * 

INPUT 
*  *  * 
*  *  * 
*  *  * 


seam# ( 1  - ' , I  1 ,  ' ) ,  regi on( 1  - ' , 12 ,  ' 
SEAM  ',11,'  IS  NOT  DEFINED.'//) 
depth) y<  =  0 . 0 ) ,  reg i on( 1  - ' , 1 2 , ' ) . 
DEPTH  must  be  <=  0.0  DEPTH  was: 
SEGMENT  RANGE  (',12, '-',12,')  IS 
BAD  OPTION: ' . 12//) 


)  .  .  . (315)  '  ) 

.  .  (F10.  1,215)'  ) 

' ,F10. 1 ) 

NOT  DEFINED . '//) 


SUBROUTINE  OFSEAM(LVL, ISEGB, ISEGE ) 

C 

C--COMPUTE  STRESSES  AND  DISPLACEMENTS  AT  OFF-SEAM  LOCATIONS. 

C 

c - 

COMMON/CONS/HW, DENSE , SR AT  10 , NSEAMS , NSEG 

COMMON/CM  1/  GAM  1 ,GAM2,Q1 ,Q2,CN1 , CN2.S1 1 , S22 , S 1 2 , C66 , B , BC66 , 

*  GM 1 2 , GM 1 2$2 , GM 1 $2 , GM2$2 , GM 1 $3 , GM2$3 , Q 1 P 1 .Q2P1 
COMMON/ CM2/  ASS (80. 5, 5) ,ASN(80,5,5) ,ANS(80,5,5) ,ANN(80,5,5) , 

1  BSS ( 80 ,5,5) , BSN ( 80 ,5,5) , BNS(80, 5 , 5 ) , BNN(80. 5 , 5 ) 

COMMON/ CM3/  DEPTH) 5) ,PS(5) ,PN(5) ,THICK(5) ,ESEAM(5) , 

*  GSEAM(5) ,MIN(80,5) 

C0MM0N/CM4/  DS ( 80 , 5 ) , DN( 80 , 5 ) 

C0MM0N/CM6/  UX ( 80 , 5 ) , UY ( 80 , 5 ) ,ESXX(80,5) , ESYY ( 80 , 5 ) , ESXY  ( 80 , 5 ) 

1  S I GXX ( 80 ,5),SIGYY( 80 ,5) ,SIGXY(80,5) ,H(5) 


Y=H( LVL ) 


o  o  o  o  o  o 


. 

, 


o  o 
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DO  500  ISEG=ISEGB, ISEGE 
UXX=0. 

UY Y=0 . 

SGY Y  =  -DENSE  *Y 
SGXX=SRATIO*SGYY 
SGX Y  =0 . 

DO  450  M  = 1 .NSEAMS 
HM  =  DE  PTH ( M ) 

DO  450  I M= 1 , NSEG 
DDS=DS ( I M , M ) 

DDN=DN ( IM,M) 

LD=2* ( ISEG-IM) 

XP  =  HW* ( LD+ 1 ) 

XM=HW*(LD-1 ) 

CALL  GFUNC( XP , HM , Y , G 1 P , G2P , G3P , G4P , G5P , G6P , G7P , G8P , G9P , 

*  G10P) 
CALL  GFUNC(XM,HM, Y , G 1 M , G2M , G3M , G4M , G5M , G6M , G7M , G8M , G9M , 

*  G10M) 
SGXX=SGXX  +  BC66*( ( G 1 P-G 1 M ) *DDS+ ( G2P-G2M ) *DDN ) 

SGY Y=SGYY  +  BC66 * ( ( G3P-G3M ) *DDS  + ( G4P -G4M ) *DDN ) 

SGXY=SGXY  +  BC66*( ( G5P-G5M ) *DDS+ ( G6P-G6M ) *DDN ) 

UXX=UXX  +  B*  ( (G7P-G7M)*DDS+( G8P-G8M ) *DDN ) 

UY Y  =UY Y  +  B*( (G9P-G9M)*DDS+( G 10P-G 10M ) *DDN ) 

450  CONTINUE 

EXX  =  S 1 1  *  SGXX  +  S 1 2*SGYY 
E YY  =  S22*SGYY  +  S12*SGXX 
EXY=0.5*SGXY/C66 
ESXX( ISEG, LVL)=EXX 
ESYY( ISEG, LVL)=EYY 
ESXY( ISEG, LVL)=EXY 
UX( ISEG, LVL)=UXX 
UY( ISEG, LVL)=UYY 
SIGXX( ISEG, LVL)=SGXX 
SIGYY( ISEG, LVL)=SGYY 
S I GX Y ( ISEG, LVL)=SGXY 
500  CONTINUE 
RETURN 
C 


END 


SUBROUTINE  ONSEAM( ISM, ISEGB , ISEGE ) 

COMPUTE  STESSES  AND  DISPLACEMENTS  AT  SEAM  LEVEL  INDICATED. 


COMMON/CONS/HW, DENSE , SR AT  10 , NSEAMS , NSEG 

COMMON/ CM2/  ASS (80, 5 , 5 ) , ASN( 80, 5 , 5 ) , ANS(80, 5 , 5) , ANN (80, 5 , 5) , 

1  BSS( 80, 5,5) , BSN ( 80 ,5,5) ,BNS(80,5,5) ,BNN(80,5,5) 

COMMON/ CM3/  DEPTH( 5 ) , PS( 5 ) , PN( 5 ) , THICK( 5 ) , ESEAM( 5 ) , 

*  GSEAM(5) ,MIN(80,5) 

COMMON/ CM4/  DS ( 80 , 5 ) , DN( 80 , 5 ) 

C0MM0N/CM5/  UXP ( 80 , 5 ) , UXM ( 80 , 5 ) , UYP ( 80 . 5 ) , UYM( 80 . 5 ) , S I GS ( 80 . 5 ) 
1  S I GN( 80 , 5 ) 


DO  20  ISEG=ISEGB , ISEGE 
UP  =  0. 

UM  =  0. 

VP  =  O. 

VM  =  0. 


C 


o  o  o  o  o  n 


oo  ooooo  no  no 


182 


SGS  =  PS(ISM) 

SGN  =  PN(ISM) 

DO  10  M  =  1 , NSEAMS 
DO  10  JM  =  1 , NSEG 
DDS=DS( UM,M) 

DDN=DN( JM , M ) 

LD  =  ISEG  -  JM 
K  =  IABS(LD)  +  1 
FAC  =  I S I GN( 1 , LD) 

SGS  =  SGS+  ASS(K,M, ISM)*DDS+FAC*ASN(K,M, ISM) *DDN 
SGN  =  SGN+F  AC  * ANS (K,M, ISM)*DDS+  ANN ( K , M , I SM ) *DDN 

TX  =  BSS(K,M, ISM)*DDS+FAC*BSN(K,M, ISM)*DDN 

TY  =  FAC*BNS(K,M, ISM)*DDS+  BNN ( K , M , ISM) *DDN 
UP  =  UP  +  TX 
UM  =  UM  +  TX 
VP  =  VP  +  TY 
VM  =  VM  +  TY 
IF  (ISM.NE.M)  GO  TO  10 
IF  ( ISEG.NE . JM)  GO  TO  10 
UP  =  UP  -  0 . 5 *DDS 
UM  =  UM  +  O . 5*DDS 
VP  =  VP  -  0.5*DDN 
VM  =  VM  +  0 . 5*DDN 
10  CONTINUE 

UXP( ISEG, I SM ) =UP 
UXM( ISEG, I SM ) =UM 
UYP ( ISEG, I SM  )  = VP 
UYM ( ISEG, I SM ) = VM 
S I GS ( ISEG, I SM ) =SGS 
SIGN( ISEG, ISM)=SGN 
20  CONTINUE 


RETURN 

END 


SUBROUTINE  GFUNC ( X , H , Y , G 1 , G2 , G3 , G4 , G5 , G6 , G7 , G8 , G9 , G 1 0 ) 

C 

--  COMPUTE  THE  TEN  G-FUNCTION  VALUES  AS  GIVEN  BY  EO . ( 37 . 33 ) - ( 37 . 45 )  C 
--  OF  CROUCH,  STEVEN  L.  C 

C 

- c 

COMMON/CM  1/  GAM  1 ,GAM2,Q1 .02.CN1 , CN2 , S 1 1 , S22 , S 1 2 , C66 , B , BC66 , 

*  GM 1 2 , GM 1 2$2 , GM 1 $2 , GM2$2 , GM 1 $3 , GM2$3 , 0 1 P 1 ,Q2P1 


YNH= Y -H 
YPH=Y+H 
YNH$2  =  YNH*  *2 
YPH$2=YPH**2 
GM 1 Y  =  GAM 1  * Y 
GM2 Y  =  GAM2  *  Y 
GM 1 X=GAM 1 *X 
GM2X  =  GAM2  *X 
GM 12X=GM12*X 
G  1  2  X  $  2  =  GM 1 2  X  *  *  2 
GM 1 X$2  =  GM 1 X  *  *2 
GM2X$2=GM2X**2 
GM 1 H=GAM 1 *H 


O  O 


' 
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GM2H  =  GAM2*Fi 
G1 YG2H=GM1 Y-GM2H 
G2 YG 1 H=GM2 Y-GM 1 H 
C 

U 1 = 1 . /(GM1X$2+YPH$2 ) 

U2= 1 . /(GM2X$2+YPH$2 ) 

GU 1 2  =  GAM 1 *U 1 -GAM2*U2 
V 1 = 1 . / ( GM 1 X$2+ YNH$2 ) 

V2  = 1  . /(GM2X$2+YNH$2 ) 

W 1 = 1 ./ (G12X$2+G2YG1H**2) 

W2= 1 . /(G12X$2+G1 YG2H**2 ) 

R1  =  -AL0G( GM 1 $2  *U 1 )/01P1 
R2=-AL0G(GM2$2*U2 )/Q2P 1 
S1  =  -ALQG( GM 1 $2  * V 1 )/Q1P1 
S2=-AL0G(GM2$2*V2)/Q2P1 
T  1  =  -AL0G( GM 1 2$2  * W 1 )/Q1P1 
T2=-AL0G( GM 12$2*W2)/Q2P1 
Z 1 =G 1 YG2H*W2 
Z2  =  G2 YG 1 H* W 1 
P=GAM2*Z2+GAM1 *  Z 1 
W 1 PW2  =  W 1 + W2 

F=CN1 *(GAM1 *V1+GAM2*V2 ) 

F  YNH  =  F  *  YNH 
C 

G1 =-YPH*(U1 /GAM 1-U2/GAM2 )-CN1 *YNH*( VI /GAM 1+V2/GAM2 ) +CN2/GM 1 2 *P 
G2=X*(GU12-F+CN2*(GM2$2*W1+GM1$2*W2) ) 

G3  =  YPH*GU 1 2  +  F  YNH-CN2  *(GAM1*Z2  +  GAM2*Z1 ) 

G4=X*(-( GM 1 $3*U 1 -GM2$3*U2 )  +  CN 1  * ( GM 1 $3* V 1 +GM2$3* V2 ) -CN2* 

1  GM 1 2$2  * W 1 PW2 ) 

G5  =  X  * ( GU 1 2  +  F -  GM 12*CN2*W1PW2) 

G6=YPH*GU12-FYNH+CN2*P 

C 

G8  =  -0. 5* ( (GAM1*R1 -GAM2*R2 )-CN1 * ( GAM  1 *S 1 +GAM2 *S2 ) +CN2 * ( T 1 +T2 ) ) 

G9  =  -0 . 5* ( (01 *R 1/GAMI -Q2  +  R2/GAM2 )+CN1 *(01 *S 1 /GAM  1 +02 * S2/GAM2 )- 
1  CN2*(Q1*T 1+Q2*T2)/( GM 12)) 

S 1 = AT  AN ( GM 1 X/ YNH ) /Q 1 P 1 
S2=ATAN(GM2X/YNH)/02P1 
T 1 =ATAN(GM12X/G2YG1H)/Q1P 1 
T2=ATAN(GM12X/G1 YG2H)/Q2P1 
R  1  =0  . 

R2  =0 . 

IF  ( YPH . EQ . O . )  GO  TO  15 
R1=ATAN(GM1X/YPH)/Q1P1 
R2=ATAN( GM2X/YPH ) /Q2P 1 
15  CONTINUE 

G7=(R1-R2)+CN1*(S1+S2)-CN2*(GAM1*T1+GAM2*T2)/( GM 1 2 ) 
G10=-(Q1*R1-Q2*R2)+CN1*(Q1*S1+Q2*S2)-CN2*(GAM2*Q1*T1+GAM1*Q2*T2)/ 

1  ( GM 1 2  ) 

C 

RETURN 

END 

***************************************************************  ******* 
SUBROUTINE  GRAPH( XX , Y Y , N2 , M , Y Y ID , XT  I TLE , NX , YT I TLE , NY . GT I TLE , I P I CT ) 


C  C 

C -  Routine  to  draw  up  to  5  curves  on  one  graph.  The  scale  C 

C -  parameters  of  the  most  deviate  curve  are  used  for  all  C 

C -  curves.  Each  curve  must  have  the  same  number  of  points.  C 

C  M  no.  of  curves  C 

C  N2  no.  of  data  points  +2  C 


o  o 


' 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


XX 

YY(  1  .  j  ) 
VYID 
XTITLE 
NX 

YTITLE 

NY 

GTITLE 


vector  of  absisca  values 

vector  of  j  '  th  ordinate  values 

vector  of  16  char,  words  describing  YY(1,j) 

vector  of  characters  to  label  x-axis. 

no.  of  chars,  in  XTITLE  (see  AXI S-*PLOTLIB ) 

vector  of  characters  to  label  y-axis. 

no.  of  chars,  in  YTITLE  (see  AXIS-*PLOTLIB ) 

vector  of  characters  used  to  label  the  graph. 

This  string  of  characters  is  terminated  by 

the  characters  "<E>". 


*****************, 


COMPLEX* 16  Y Y I D ( 1 ) 

DIMENSION  XX (  1  )  , YY(N2,  1  )  , ISYMBL(5) 

INTEGER  XTITLE (  1  )  , YTITLE( 1 ) ,GTITLE( 1 ) 

DATA  I SYMBL/ ' a<E> ' , 'b<E>' , ' c<E> ' , 'd<E>' , ' e<E> ' / 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

******  *Q 


N=N2-2 

CALL  I GBGNS ( I  PICT  .  ' WIND ' .-0.5,  11  . 5 , -0. 5 .8 . 5  ) 
CALL  SCALE(XX,9.0,N,  1  ) 


find  common  scaling  parameters  for  all  curves 
MN2=M*N2 

CALL  SCALE(YY( 1 , 1 ) .8.0.MN2, 1 ) 

AMIN= YY ( 1.M+1) 

ADELT=YY ( 2 , M+ 1 ) 


set  scale  parameters  for  each  curve 
DO  50  d= 1 ,M 

Y Y ( N+ 1 ,U)=AMIN 
YY(N+2,U)=ADELT 
50  CONTINUE 


draw  the  curves. 

ICHAR=0 
DO  150  U= 1 , M 

CALL  LINE(XX,YY( 1 , J) ,N,  1  ,2, I CHAR) 

I CHAR  =  I CHAR+ 1 
150  CONTINUE 

plot  first  the  X-axis  followed  by  the  Y-axis. 

CALL  AXIS2(0. .0. .XTITLE , -NX , 9 . , 0 . , XX ( N+ 1 ) ,XX(N+2 ) . 1 . ) 

CALL  AX  I S ( 0 .  ,0.  .YTITLE, NY, 8.  ,90.  , AMIN , ADELT ) 

describe  the  curves. 

CALL  I GMA ( 9 . 5 , 8 . 0 ) 

DO  200  U= 1 , M 

CALL  I GTXT ( '<ASCL>' , . 18, '<F0NT>' , 'PLOT' , ISYMBL(d) ) 

CALL  I GTXT ( ' <FONT> ' , ' 7ASC ' , '  =  <E> ' , Y Y I D ( J ) , ' <CRLF><CRLF><E > ' ) 
200  CONTINUE 

label  the  graph. 

CALL  IGMA( 1 .5,8.0) 

CALL  IGTXT( '<ASCL>' ,. 18, GTITLE) 

CALL  I GENDS ( I P I CT ) 

RETURN 
END 

£****************************************************★*************0 
C  C 


non  n  n 
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SUBROUTINE  GETYY(M,Y1 , IDY 1 ,Y2, IDY2,Y3, IDY3,N2,YY,YYID) 

Routine  to  move  the  curve  data  into  columns  of  a  single 
matrix  for  the  graph  plotting  routine,  GRAPH.  Curve  label¬ 
ling  text  is  also  assembled  for  each  curve. 

M  number  of  curves  to  be  assembled. 

Yi  vector  containing  i'th  curve  ordinate  values 

IDYi  16  char,  word  describing  i'th  curve 
N2  no.  of  data  points  +  2 

Y Y (  1  , j  )  vector  of  j'th  ordinate  values  (output) 

YYID  vector  of  16  char,  words  defining  curves ( output ) 


******************* 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

*c 


COMPLEX* 16  IDY 1 , IDY2, IDY3, YYID( 1 ) 
DIMENSION  Y1(1),Y2(1),Y3(1),YY(N2,1) 


N=N2-2 
Y  Y  I D  (  1  )  =  I D  Y  1 
DO  10  1  =  1, N 

YY( I ,  1  )  =  Y1  (  I  ) 

10  CONTINUE 

IF  (M.LT.2)  GOTO  40 
YYID(2)=IDY2 
DO  20  1=1 ,N 

Y Y (  I  ,  2  )  = Y2 ( I ) 

20  CONTINUE 

IF  (M.LT.3)  GOTO  40 
YYID(3)=IDY3 
DO  30  1= 1 ,N 

Y  Y  (  I  ,  3  )  =  Y  3  (  I  ) 

30  CONTINUE 


40  DO  50  >J=  1  ,M 

Y Y ( N+ 1 . J)=YY(N, J) 
YY(N+2,U)=YY(N, J) 
50  CONTINUE 
C 


RETURN 

END 

C= ======================================================= ===C 

SUBROUTINE  ME SSGE ( I STRNG , X , Y  ) 

C  C 

C  —  to  place  a  message  " I STRNG " ,  starting  at  ABS(X.Y)  C 

C  C 

c= ========================================================= =c 


DIMENSION  I STRNG(  1  ) 

CALL  IGBGNS( 'MESSAGE '  ) 

CALL  IGMA(X.Y) 

CALL  I GTXTH ( I STRNG ) 

CALL  IGENDS( 'MESSAGE ' ) 

RETURN 

END 

C 

INTEGER  FUNCTION  I TEM ( I CMD , X , Y ) 

C 

-  Function  to  place  menu  Item  at  position  X,Y  and  C 

-  return  a  picture  name  for  the  item.  C 


, 


o  o  o  o 
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c 

DIMENSION  I CMD ( 1 ) 

I T  EM= I GBGNS ( O ) 

CALL  IGMA(X.Y) 

CALL  IGTXT(ICMD) 

CALL  I GENDS (ITEM) 

RETURN 

END 


C 


C 


C 


C 


SUBROUTINE  STATS(IPICT) 

C 

Routine  to  plot  the  statistics  describing  the  c 

multiple  seam  mine.  q 

C 

COMMON/CONS/HW, DENSE , SRATIO , NSEAMS , NSEG 

COMMON/ CMO/  ICARD(20),EXX,EYY,EZZ, VYX , VZX , VYZ , GXY 

COMMON/ CM3/  DEPTH( 5 ) , PS  (  5  ) , PN( 5 ) ,THICK(5) , ESEAM(5) , 

*  GSEAM(5) ,MIN(80,5) 

CALL  I GBGNS ( I P I CT , ' WIND ' ,0. , 1 . ,0. , 1 . ) 

CALL  I GMA ( 1 .0,0.0) 

CALL  IGDR(- .998,0.0) 

CALL  I GDR (0.0, 1.0) 

CALL  I GMA (0.01 , .94) 

CALL  IGFMT(ICARD, 'A' ,35) 

CALL  I GMA (0.01 , .84 ) 

CALL  I GT  XT ( ' N<BSUB>SEAM<ESUB>  =  <E> ' ) 

CALL  IGFMT(NSEAMS, ' I ' ,3) 

CALL  I GMA (0.01 , .78) 

CALL  IGTXT( 'N<BSUB>SEG  <ESUB>=<E>') 

CALL  I GFMT ( NSEG , '  I  '  ,  3) 

CALL  I GMA (0.01 ,  .12) 

CALL  I GTXT ( ' WD<BSUB>HALF<ESUB>=<E> ' ) 

CALL  I GFMT ( HW , 'F' ,6,2) 

CALL  I GMA (0.01 , . 66 ) 

CALL  I GTXT (  ' WT<BSUB>UNI T<E SUB>  =  <E> '  ) 

CALL  IGFMT(DENSE, 'F' ,7,3) 

CALL  I GMA (0.01 , .6) 

CALL  I GTXT ( ' S I G<BSUB>XX<E SUB>/S I G<BSUB>Y Y<ESUB>=<E> ' ) 
CALL  IGFMT(SRATIO, 'F' ,7,3) 

CALL  I GMA (  . 4 ,  . 84  ) 

CALL  I GTXT ( ' E<BSUB>XX<ESUB>  =  <E> '  ) 

CALL  I GFMT ( EXX , 'E' ,7,2) 

CALL  I GMA (  . 4 ,  . 78  ) 

CALL  I GTXT ( ' E<BSUB>YY<ESUB>  =  <E> '  ) 

CALL  I GFMT ( E Y Y , 'E' ,7,2) 

CALL  I GMA ( .4, .72) 

CALL  I GTXT ( ' E<BSUB>ZZ<ESUB>  =  <E>  '  ) 

CALL  IGFMT(EZZ, 'E' ,7,2) 

CALL  I GMA ( .4, .66) 

CALL  I GTXT (  ' G<BSUB>XY<ESUB>  =  <E> '  ) 

CALL  I GFMT (GXY, 'E' ,7,2) 

CALL  I GMA ( .75, .84) 

CALL  I GTXT (  ' V<BSUB>XY<ESUB>  =  <E>  '  ) 


o  o 


. 
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CALL  I GFMT ( VYX ,  'F  '  ,3,2) 

CALL  I GMA (  .75,  .  78  ) 

CALL  I GTXT ( ' V<BSUB>ZX<ESUB>=<E> '  ) 
CALL  IGFMT(VZX, 'F ' ,3,2) 

CALL  I GMA (  . 75 ,  . 72  ) 

CALL  I GTXT ( ' V<BSUB>YZ<ESUB>  =  <E> '  ) 
CALL  IGFMT ( VYZ , 'F' ,3,2) 


CALL  I GMA (0.01 , .5) 

CALL  I GTXT ( ' SM#  E<BSUB>SEAM<ESUB><E> '  ) 

CALL  I GTXT (  G<BSUB>SE AM<ESUB>  THICK  DEPTH<E> ' ) 

CALL  I GMA (0.01 , .46) 

DO  20  U  = 1 .NSEAMS 


CALL 

I GTXT ( ' <CRLF><E> ' ) 

CALL 

I GFMT ( d , ' I ' ,2) 

CALL 

I GTXT (  '  <E > '  ) 

CALL 

IGFMT(ESEAM(U) , 'E' ,8,3) 

CALL 

I GTXT ( '  <E> '  ) 

CALL 

I GFMT ( GSE AM ( J ) , 'E' ,8,3) 

CALL 

I GTXT (  '  <E>  '  ) 

CALL 

IGFMT(THICK(J) , ' F ' , 5 ) 

CALL 

I GTXT ( '  <E> ' ) 

CALL 

I GFMT ( DEPTH ( J ) , 'F'  ,6,  1  ) 

CONTINUE 

CALL  IGENDS( IPICT ) 

RETURN 

END 

SUBROUTINE  MINEX(IMINE) 

C 

Routine  to  plot  a  cross-section  of  the  multiple  seam  mine  C 

that  is  under  investigation.  The  same  x-axis  as  used  for  c 

the  graphs  is  used  here  to  ease  interpretation.  Also  the  C 

seam  segmentation  is  drawn  as  an  additional  aid.  c 


COMMON/CONS/HW, DENSE , SR AT  1 0 , NSE AMS , NSEG 
COMMON/ CM3/  DEPTH( 5 ) , PS ( 5 ) , PN( 5 ) . THICK) 5 ) . ESEAM( 5 ) , 
*  GSEAM(5) ,MIN(80,5) 

DIMENSION  X ( 4  )  , Y(4) 


I M I NE  =  I GBGNS ( 0 ,  'WIND'  ,-.5, 9. 5, -.5, 9. 5) 
Y ( 1  )  =  -VMIN( DEPTH, NSE AMS) 

Y( 2 )=-VMAX(DE PTH, NSEAMS) 


verify  that  the  scaling  in  the  y-direction 
seams  appear  about  1/20'th  of  the  vertical 
AVHO=VAVG( THICK, NSEAMS) 

RANGE=ABS( Y ( 1 )-Y(2) ) 

IF  (AVHO. LE . RANGE/20. )  GOTO  10 
HRNGE  =  ( Y ( 1 ) +Y ( 2 ) )* .5 
Y ( 1 ) =HRNGE+ 10 . * AVHO 
Y ( 2 ) =HRNGE -  1 0 . *AVHO 
10  CALL  SCALE ( Y , 9 . 0 , 2 , 1) 

YMIN  =  Y ( 3  ) 

YDELT  =  Y ( 4 ) 


i s  such 
sea  1 e . 
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use  the  same  scaling  in  the  x-direction  as  used  for 
X( 1 ) =HW*NSEG 
X(2)  =  -X(  1  ) 

CALL  SCALE ( X , 9.0,2, 1 ) 

XMIN  =  X ( 3  ) 

XDELT  =  X ( 4 ) 


the  graphs. 


draw  a  segmentation  axis,  annotating  every  5'th  segment  along 
W  =  2 . *HW 
NSEG  1 =NSEG+ 1 
DT I C  =  5 . * W/XDE  LT 
AX=(X(2) -HW-XMIN )/XDELT 
AXLEN=(W+NSEG1 )/XDELT 
AXDELT=NSEG1/AXLEN 
CALL  AX2EP(DTIC,2,0,4) 

CALL  AXIS2(AX,9. , 'SEAM  SEGMENTATION' 
draw  X  and  V  axes . 

CALL  AX2EP (  1  .  ,2,0,0,D,D,D,D,D,D,D,D) 

CALL  AXIS2(0. ,0. , 'HORIZONTAL  D I  STANCE '  , -  1 9 , 9 . 0 , 0 . 

CALL  AX  I S2 ( O .  ,0.  ,  'DEPTH'  ,5,9.0,90.  , YMI N , YDE  LT ,  1  .  ) 


seam 


1  7  ,  AXLEN , 0 .  .0.  , AXDELT , DT I C ) 


, XMIN , XDELT , 1 . ) 


draw  each  seam. 

DO  195  d=  1  , NSE AMS 
DE  L Y  =  THI CK (  d  ) 

HE IGHT  =DELY/YDELT 
AX  =  X ( 2  ) 

VY=-(DEPTH(d)+DELY* . 5+ YMIN ) /YDE LT 
I  MI N  =  MIN(  1  ,  d) 

ICNT  =  2 

DO  175  1=2, NSEG 

IF(MIN( I , J) . EQ. IMIN)  GOTO  150 
DE  LX  =  I CNT  *HW 

130  CALL  IBAR( ( AX-XMIN )/XDELT , VY .DELX/XDELT .HEIGHT , IMIN, 10) 

AX=AX+DELX 
IMIN=MIN( I , d) 

I CNT  =  2 
GOTO  175 

150  I CNT  = I CNT  +  2 

175  CONTINUE 

DE  LX  =  I CNT  *HW 

CALL  I BAR ((AX-XMIN)/XDELT,VY, DELX/XDELT, HE IGHT. IMIN, 10) 
label  the  seam  by  number. 

CALL  I GMA ( (AX+DELX-XMIN) /XDELT, VY) 

CALL  IGTXT( '<ASCL>' .HEIGHT* .65, '  SM#<E>') 

CALL  I GFMT ( d ,  ' I '  ) 

195  CONTINUE 

CALL  IGENDS( IMINE ) 


RETURN 

END 

********************************************************************** 


SUBROUTINE  I BAR ( VX , VY , DE LX , DELY , ITYPE.LPB) 

C  C 

C -  Routine  to  draw  a  seam  element  as  a  mined  panel  or  as  an  C 

C -  unmined  pillar  at  point  (VX.VY).  C 

C  VX.VY  lower  left  corner  of  element  C 

C  DELX.DELY  element  width  and  height  C 


o  o 
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c 

c 

c 

c 

Q********* 

c 


I  TYPE  =  1 
I  TYPE  =  2 
LPB 

******** 


mined  e 
unm i ned 
no .  of 

******** 


lenient,  no  hatching 
element,  hatching  1 
hatch  lines  per  box 

******************** 


s  performed 


*********************** 


C 

C 

c 

c 

*c 


CALL  IGMA(VX.VY) 

CALL  IGDR(0.0,DELY ) 

CALL  IGDR(DELX ,0.0) 

CALL  IGDR(0.0, -DELY) 

CALL  IGDR(-DELX.O.O) 

IF  (ITYPE.EO. 1)  RETURN 
DH  =  DE  L Y/ LPB 
DX=DELX 
DY  =0 . 0 

DO  10  J=2 , LPB 

CALL  I GMR (0.0, DH ) 

CALL  IGDR(DX.O.O) 
DX=-DX 


10  CONTINUE 
C 


RETURN 

END 

REAL  FUNCTION  VAVG(VEC.N) 

DIMENSION  VEC ( 1  ) 

VAVG= VEC ( 1 ) 

IF  (N. LE . 1 )  RETURN 
DO  10  J=2 , N 

VAVG=VAVG+VEC(U) 

10  CONTINUE 

VA VG= VA VG/N 
RETURN 
END 

REAL  FUNCTION  VMIN(VEC.N) 

DIMENSION  VEC ( 1 ) 

VMIN  =  VEC(  1  ) 

IF  (N.LE.1)  RETURN 
DO  10  U=2,N 

IF  (VEC(U) .LT.VMIN)  VMIN=VEC(J) 
10  CONTINUE 
RETURN 
END 

REAL  FUNCTION  VMAX(VEC.N) 

DIMENSION  VEC ( 1  ) 

VMAX  =  VEC (  1  ) 

IF  (N.LE.1)  RETURN 
DO  10  J=2 , N 

IF  (VEC(d) .GT.VMAX)  VMAX=VEC(U) 
10  CONTINUE 
RETURN 
END 


■ 
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Program  DDSEAMS  sample  output 


RUN  02  TEST  CRACK  NSEG=40 
NUMBER  OF  SEPARATE  SEAMS  =  1 

NUMBER  OF  SEGMENTSIN  EACH  SEAM  =  40 

HALF-WIDTH  OF  SEGMENTS  =  2.00 

UNIT  WEIGHT  OF  ROCK  =200.00 


RATIO  OF  PRE-EXISTING  STRESSES  SIGXX/SIGYY  =  0.0 


ELASTIC 

CONSTANTS 

EXX  = 

0. 7000E+08 

EYY  = 

0. 7000E+08 

EZZ  = 

0. 70O0E+08 

GXY  = 

0. 3900E+07 

VYX  = 

0.25 

VZX  = 

0.25 

VYZ  = 

0.25 

ELASTIC 

PARAMETERS 

GAM  1 
GAM2 
01 
02 


1 .00 
9 . 24 
4 . 29 
0.23 
48 . 39 
0.02 


SPECIFICATIONS  FOR  SEAM  NUMBER  1 

SEAM  THICKNESS  =  10.00 

DEPTH  =  2000.00 

MODULUS  OF  ELASTICITY  =  0.700E+08 
SHEAR  MODULUS  =  0.390E+07 

MINING  PATTERN  (1  =  MINED,  2  =  UNMINED) 

2222222222111111111111111111112222222222 


SOLUTION  DATA  FOR  GAUSS-SIEDEL  ITERATION  PROCESS 


NUMBER  OF  ITERATIONS:  31 
RELAXATION  FACTOR:  1.60 
ATTAINED  TOLERANCE:  0.659E-05 


' 

STRESSES,  DISPLACEMENTS  AND  DISPLACEMENT  DISCONTINUITY  COMPONENTS  AT  SEAM  H 1 
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APPENDIX  2 


CRACK  PROGRAM 
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Program  CRACK  listing 


READ  (4.100)  W , ROE , PZ , P Y , PX 

100  FORMAT  (5F10.0) 

WRITE  (7,201)  W , ROE , PX , P Y , PZ 

201  FORMAT  (////// 1 X CRACK  SOLUTION  PROGRAM  CRACK . S  OUTPUT'// 

+10X, 'WIDTH  =  ' ,F10. 1/ 

+10X, 'POISSONS  RATIO  =  '.F10.4/ 

+5X , ' PX  =' .E10.3.5X, 'PY  =  ',E10.3,'PZ  =  '.E10.3/// 

+  5X,  'X'  , 12X,  'SIGX'  ,  1  IX,  'SIGY'  ,  1 IX,  'SIGZ'/) 

1  READ  (4,101)  X 

101  FORMAT  (F10.0) 

IF  (X.EO. 999.0)  GO  TO  998 

XLAM  =  1 .0/SQRT( 1 .0~W*W/(4 .0*X*X) ) 

SIGX  =  (XLAM-1 .0)*PZ  +PX 

SIGY  =  PY  +  2.0*R0E*(XLAM-1 .0)*PZ 

SIGZ  =  XLAM*PZ 

WRITE  (7,202)  X , S I GX , S I GY , S I GZ 

202  FORMAT  ( 4 ( E 10 . 3 . 5X ) / ) 

GO  TO  1 

998  WRITE  (7,203) 

203  FORMAT  (///// 

+5X, 'Z' , 12X, 'SIGX' , 1 IX, 'SIGY' , 1 IX, 'SIGZ'/) 

2  READ  (4, 102)  Z 

102  FORMAT  (F10.0) 

IF  (Z.EO. 999.0)  GO  TO  999 

ZLAM  =  1 .0/SQRT( 1 ,0+W*W/(4 .0*Z*Z) ) 

SIGX  =  (2 ,0*ZLAM-ZLAM*ZLAM*ZLAM-1 .0)*PZ+PX 

SIGY  =  PY+2.0*ROE*(ZLAM-1  .0)*PZ 

SIGZ  =  ZLAM*ZLAM*ZLAM*PZ 

WRITE  (7,202)  Z , S I GX , S I GY , S I GZ 

GO  TO  2 

999  STOP 
END 


' 
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Program  CRACK  sample  output 


C-7ACK  SOLUTION  PROGRAM  CRACK  .  S  OUTPUT 

WIDTH  =  80.0 

POISSONS  RATIO  =  0.2500 

Px  =  0.0  PY  =  0.0  PZ  =  0 . 400E+06 


X 

SIGX 

SIGY 

SIGZ 

0. 4 10E+02 

0. 142E+07 

0.71 1 E  +06 

0. 182E+07 

0. 430E+02 

0. 690E+06 

0. 345E+06 

0. 109E+07 

0 . 450E+02 

0. 473E+06 

0. 237E+06 

0 . 873E+06 

0. 470E+02 

0. 3G2E+06 

0. 18 1E+06 

0. 762E+06 

0. 490E+02 

0. 293E+06 

0. 146E+06 

0. 693E+06 

0.510E+02 

0. 245E+06 

0. 122E+06 

0 . 645E+06 

0. 530E+02 

0. 2 1 0E  +06 

0. 105E+06 

0 . 6 10E+06 

0. 550E+02 

0. 183E+06 

0.91 4E+05 

0. 583E+06 

0. 570E+02 

0. 161E+06 

0. 807E+05 

0.561E+06 

0. 590E+02 

0 . 144E+06 

0. 721E+05 

0.544E+06 

0 . 6 10E+02 

0. 130E+06 

0 . 649E+05 

0. 530E+06 

0.630E+02 

0. 1 18E+06 

0 . 589E+05 

0. 5 18E+06 

0.650E+02 

0. 107E+06 

0. 537E+05 

0. 507E+06 

0. 670E+02 

0. 986E+05 

0. 493E+05 

0 . 499E+06 

0. 690E+02 

0. 909E+05 

0. 455E+05 

0. 491E+06 

0. 7 10E+02 

0. 841E+05 

0. 42 1E+05 

0.484E+06 

0. 730E+02 

0. 782E+05 

0. 391E+05 

0. 478E+06 

0. 750E+02 

0. 729E+05 

0. 364E+05 

0 . 473E+06 

0. 770E+02 

0.681E+05 

0.341E+05 

0. 468E+06 

0. 790E+02 

0 . 639E+05 

0. 319E+05 

0. 464E+06 

Z 

SIGX 

SIGY 

SIGZ 

0. 200E+01 

-0. 360E+06 

-0. 190E+06 

0. 498E+02 

0. 500E+01 

-0. 302E+06 

-0. 175E+06 

0. 763E+03 

0. 100E+02 

-0. 212E+06 

-0. 151E+06 

0.571E+04 

0. 250E+02 

-0. 356E+05 

-0. 940E+05 

0. 596E+05 

0. 350E+02 

0. 126E+05 

-0.683E+05 

0. 1 14E+06 

0 . 450E+02 

0. 309E+05 

-0.505E+05 

0. 167E+06 

0. 650E+02 

0. 342E+05 

-0. 297E+05 

0. 247E+06 

0 . 950E+02 

0. 242E+05 

-0. 157E+05 

0.313E+06 

. 
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APPENDIX  3 


PSTRESS  PROGRAM 


o  n 
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Program  PSTRESS  listing 


DIMENSION  C(5) ,PHI(5),RPHI(5) 

C ( 1 )  =  2.3 
PHI ( 1 )=  40.0 
C ( 2 )  =  0.55 
PHI ( 2 )  =  40.0 
C ( 3  )  =  2.3 
PHI ( 3  )  =  33.0 
C ( 4 )  =  0.0 
PHI ( 4  )  =  0.0 
C ( 5 )  =  0.0 
PHI ( 5 )  =  0.0 
BOT 1  =  -148.75 
BOT  2  =  -151.25 
BOT 3  =  -1000.0 
BOT  4  =  -9999.0 
IPAGE  =  0 
READ  (4 , 1001 )  N 
WRITE  (7,2001) 

DO  50  J  = 1 , 5 

50  RPHI(d)  =  PHI(J)*3. 141593/180.0 
DO  100  1=1 ,N 

READ (4 , 1000)  NUMEL .GRIDY ,GRIDX,DX,DY,EX,EY,EXY,X,Y,XY 
CALL  STRESS(X,Y,XY,PM1 , PM2 , TV3 , THETA , I  FLAG) 

IF  (GRIDY .GE . B0T1 )  MAT  = 1 

IF  (GRIDY . LE . BOT 1 . AND . GRIDY . GE . B0T2 )  MAT  =  2 
IF  (GRIDY . LE .B0T2 . AND . GRIDY .GE .B0T3)  MAT=3 
IF  (GRIDY. LE.B0T3. AND. GRIDY. GE.B0T4)  MAT  =  4 
IF  (GRIDY. LE.B0T4)  MAT  =  5 
IF  ( I  FLAG . EQ . O )  GO  TO  20 
SC  =  1.0 
GO  TO  19 

20  CALL  STABIL(PM1 , PM2 , RPHI ( MAT ) , C ( MAT ) , SC ) 

19  WRITE  (8,2000)  GR I DY , GR I DX , DX , D Y , X , Y , X Y , PM  1 , PM2 , THETA , TV3  SC 
+  C(MAT),PHI( MAT ) 

IPAGE  =  IPAGE  +  1 
IF  (IPAGE. NE. 58)  GOTO  3 
IPAGE  =  1 
WRITE  (7,2001) 

3  WRITE  (7,2002)  GR I DY , GR I DX , DX , DY , X , Y , XY , PM  1 , PM2 , THETA , TV3  SC 
+  MAT 

100  CONTINUE 

1000  FORMAT  (I5.10F1O.0) 

1001  FORMAT  (15) 

2000  FORMAT  ( 2F 10 .  1 , 1 2 ( 1  PE  10 . 3 ) ) 

2001  FORMAT  ( 1H1/7X, 1HY.9X, 1 HX ,4X,6HX-DISP,4X,6HY-DISP,4X,6HSIG  XX. 4X 
+  6HSIG  Y Y , 4X , 6HT AU  XY , 5X , 4HSIG1 , 6X , 4HSIG2 , 6X , 5HTHETA , 4X , 

+  6HTAUMAX , 4X , 7HSTABLTY , 3X , 3HMAT/ ) 

2002  FORMAT  ( 2F 10 .  1 ,  10( 1  PE  10 . 2 ) , I  5 ) 

STOP 
END 

SUBROUTINE  STRE SS ( X , Y . X Y . PM  1 , PM2 , TV3 , THETA , I F LAG ) 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  c 

C  SUBROUTINE  STRESS  -  DETERMINES  THE  PRINCIPAL  STRESSES  AND  C 

C  THEIR  DIRECTIONS  C 

C  C 

C  INPUT  c 

C  C 


INPUT 


. 


WO)  -1.(0  0  0  0 
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C  X  =  SIGMA-X  STRESS 

C  V  =  SIGMA-Y  STRESS 

C  XY  =  SHEAR  STRESS 

C 

C  OUTPUT 


C  PM  1  =  MAJOR  PRINCIPAL  STRESS  C 
C  PM2  =  MINOR  PRINCIPAL  STRESS  C 
C  TV3  =  MAXIMUM  SHEAR  STRESS  C 
C  THETA  =  DIRECTION  OF  MAJOR  PRINCIPAL  STRESS  C 
c  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

C - CALCULATE  PRINCIPAL  STRESSES 

C 

I  FLAG  =  O 

IF  (X .NE .0.0. OR . Y .NE .0.0. OR .XY .NE .0.0)  GO  TO  10 

PM  1  =  0.0 

PM2  =  0.0 

THETA  =  0.0 

TV3  =  0.0 

I  FLAG  =  1 

GO  TO  9 

10  TV  1 =0 . 5* ( X+Y ) 

TV2=0 . 5* ( X-Y ) 

TV3=SQRT(TV2**2+XY**2) 

PM  1 =TV 1 +TV3 
PM2  =  TV 1 -TV3 

CALCULATE  DIRECTION  OF  PRINCIPAL  STRESSES 

I F ( T V3 )  1,2,1 
THETA=0. 785398 
GO  TO  3 
TV  1 =TV2/TV3 

I F ( ( 1  . O- ABS ( TV  1 ) ) . GT .0.0001 )GOTO  6 
I F ( T V 1 . GT .0.0)  THETA=0 . 0 
IF(TV1 . LT .0.0)  THETA= 1 . 570796 
GO  TO  3 

THETA=ARC0S(TV1 )*0. 5 
IF(XY . LT .0.0)  THET  A  =  -THETA 
TV  1 =PM 1 

IF (ABS(PM1 ) .GE . ABS(PM2) )GOTO  5 
PM  1 =PM2 
PM2  =  TV 1 

THETA=THETA+1 . 570796 
5  THETA  =  THETA*(  180. /3.  141593) 

9  CONTINUE 
RETURN 
END 

SUBROUTINE  STABIL  ( PM  1 , PM2 , PHI , C . SC ) 

PTOT  =  PM1+PM2 
PD IFF  =  PM1-PM2 

SC  =  1 .0-PDIFF/(SIN(PHI )*(PT0T+2 .0*TAN(PHI ) ) ) 

RETURN 

END 


DISP  Y-DISP  SIG  XX  SIG  YY  TAU  XY  SIG1  SIG2  THETA  TAUMAX  STABLTY  MAT 
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GEOTECHNICAL  TESTING  PROCEDURES 
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UNIAXIAL  COMPRESSION 

This  test  was  undertaken  on  specimens  of  full  core  diameter 
(2.6  inch).  The  specimens  were  cut  with  a  diamond  saw,  then 
the  ends  were  lapped,  using  a  lapping  table  with  #80  and 
#220  grit.  Specimen  width  to  height  ratios  between  2.0  and 
2.5  were  prepared  (except  for  one  mudstone  specimen  with  a 
ratio  of  1.6).  The  ends  were  lapped  so  that  they  were  plane 
and  parallel  to  within  0.002  inch. They  were  then  air  dried 
and  electrical  resistance  strain  gauges  were  attached  to 
each  specimen  using  epoxy  cement.  Two  pairs  were  placed 
around  the  central  circumference,  one  pair  for  axial  strain, 
and  the  other  for  ci rcumferent i a  1  strain. 

The  specimen  was  placed  in  an  MTS  servo  controlled 
stiff  testing  machine  (Plate  2),  with  a  spherical  seat  on 
the  bottom  and  the  machine  spherical  seat  on  the  top.  Strain 
indicators  were  connected  to  the  matched  pairs  of  gauges  in 
a  full  bridge  configuration.  The  specimen  was  loaded  at  a 
constant  machine  stroke  rate,  producing  specimen  failures 
within  approcimately  10  minutes.  Strain  readings  were  taken 
at  load  increments  of  1000  or  2000  pounds,  depending  upon 
the  expected  failure  load. 

Axial  and  lateral  strain  were  plotted  against  stress 
(Appendix  7)  and  the  elastic  moduli  were  estimated  from 
these  graphs.  The  Young's  Modulus  was  estimated,  by 
measuring  the  slope  of  the  axial  stress/strain  curve,  at  the 
stress  value  equal  to  50%  of  the  failure  stress.  The 
circumferential  strain  is  negative  (geotechnical  sign 
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convention)  and  therefore  the  sign  has  been  reversed  to 
enable  both  strains  to  be  plotted  on  the  same  graph.  The 
c i rcumf erent i a  1  modulus  is  estimated  by  measuring  the  slope 
of  the  circumferential  stress/strai n  curve  at  50%  of  the 
failure  stress.  Poisson's  Ratio  is  then  calculated  from 
(International  Society  of  Rock  Mechanics,  Commission  on 
Standardisation  of  Laboratory  and  Field  Tests,  1979): 

Poisson's  Ratio  =  Young's  Modulus/Circumferential 

Modulus 


TRIAXIAL  COMPRESSION 

Specimens  for  this  test  were  undercored  to  1.6  inch  diameter 
using  a  bench  mounted  diamond  drill  with  air  or  water  flush 
(depending  upon  rock).  The  undercored  specimens  were  cut 
with  a  diamond  saw  to  produce  width  to  height  ratios  between 
2.0  and  2.5.  This  was  not  possible  for  some  mudstone  and 
coal  samples,  as  they  did  not  withstand  the  undercoring 
process  very  well.  The  ends  were  lapped  to  be  plane  and 
parallel  to  0.002  inch.  The  j_n  situ  moisture  content  was 
preserved  as  much  as  possible  by  sealing  samples  in  plastic 
bags  prior  to  testing. 

Immediately  prior  to  testing  the  specimens  were 
unsealed  and  placed  in  the  triaxial  cell.  The  Hoek- Frank  1 i n 
triaxial  cell  consists  of  a  metal  jacket  (capable  of 
withstanding  internal  hydraulic  pressures  of  10,000  psi ) ,  a 
flexible  sleeve  (to  provide  a  seal  and  protect  the  specimen 
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from  the  hydraulic  fluid),  and  spherical  seats  at  both  ends. 
The  triaxial  cell  pressure  is  provided  by  an  intensifier 
unit  that  is  pressurised  from  a  compressed  nitrogen 
cy 1 i nder . 

The  cell  was  placed  in  the  MTS  servo  controlled  stiff 
testing  machine  (Plate  2)  and  held  in  place  by  a  small  axial 
load.  The  cell  was  then  pressurised  to  the  confining 
pressure  required  for  the  test.  For  the  higher  confining 
pressures  this  was  done  in  stages,  in  order  to  avoid 
specimen  failure  as  a  result  of  the  confining  pressure.  At 
each  stage  the  axial  load  was  increased  to  produce 
hydrostatic  stress  conditions.  The  specimen  was  then  loaded 
axially  at  a  constant  machine  stroke  rate,  producing 
specimen  failures  within  approximately  10  minutes. 

After  failure  the  specimen  was  removed  and  weighed, 
then  it  was  oven  dried  at  105°  Celsius  to  calculate  the 
moisture  content. 

BRAZILIAN  DISC 

Disc  specimens  were  prepared  by  cutting  cylindrical  samples 
perpendicular  to  the  long  axis  of  the  core.  Iri  situ  moisture 
content  was  preserved  as  much  as  possible  by  sealing 
specimens  in  plastic  bags  prior  to  testing. 

The  specimens  were  tested  in  a  5000  Kg  stepless 
compression  test  machine,  with  flat  platens  and  a  spherical 
seat  at  the  bottom.  Loading  was  carried  out  along  a 
diameter,  and  the  load  at  failure  was  estimated  from  the 
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dial  gauge  reading  on  a  proving  ring  incorporated  into  the 
machine.  After  failure  the  specimen  was  weighed,  oven  dried 
at  105°  Celsius,  and  re-weighed  to  calculate  the  moisture 
content . 

The  tensile  strength  was  obtained  from  the  following 
equation  (International  Society  for  Rock  Mechanics, 
Commission  on  Standardisation  of  Laboratory  and  Field  Tests, 
1978)  : 


SigT  =  0 . 636 ( P/D ) 

SigT  =  Tensile  strength  (MN/m2) 
P  =  Maximum  load  (N) 

D  =  Diameter  of  specimen  (mm) 


DIRECT  SHEAR 

Direct  shear  specimens,  approximately  2cm  thick,  were 
prepared  by  cutting  parallel  to  the  discontinuity  plane  to 
be  tested.  During  the  preparation  process  the  j_n  situ 
moisture  content  was  preserved  as  much  as  possible  by 
covering  the  specimens  with  a  moist  cloth.  The  specimens 
were  mounted  into  each  half  of  the  shear  box,  using  a 
sulfaset  anchorage  mix  that  expands  slightly  upon  setting. 
The  discontinuity  plane  was  arranged  so  that  it  was  parallel 
to  the  plane  of  shear. 

The  shear  box  was  mounted  into  a  Farnell  direct  shear 
apparatus  that  could  provide  a  constant  shear  rate.  Normal 
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load  was  applied  by  means  of  a  suspended  loading  platform 
and  assorted  weights.  Normal  and  shear  displacements  were 
measured  by  linear  variable  transducers,  and  shear  load  by 
strain  gauges  attached  to  a  proving  ring.  All  measurements 
were  continuously  recorded  on  an  X-Y  plotter. 

The  specimen  was  sheared  approximately  5mm  with  a 
normal  load  of  300  N,  then  returned  to  its  intial  position 
with  no  normal  load.  This  procedure  was  repeated  with  normal 
loads  of  1350  N,  2850  N  and  4350  N.  Peak  and  residual 
friction  angles  were  calculated  from  the  shear  and  normal 
loads  recorded,  with  adjustment  made  for  the  measured 
di 1 atancy . 

SLAKE  DURABILITY 

Slake  durability  tests  were  carried  out  using  the  standard 
apparatus  and  procedure  (Gyenge  &  Herget ,  1977).  The  sample 
was  divided  into  ten  pieces  with  rounded  corners,  each  of 
approximately  50g  weight.  These  were  oven  dried  at  105° 
Celsius  and  weighed.  They  were  then  placed  in  the  standard 
test  drum  and  mounted  in  the  testing  apparatus.  The  drum  was 
partially  submerged  in  tap  water  at  room  temperature,  and 
rotated  at  20  rpm  for  10  minutes.  The  sample  was  removed, 
and  oven  dried  before  being  re-weighed.  This  procedure  was 
repeated  and  the  final  weight  after  the  second  cycle  was 
obta i ned . 

The  slake  durability  index  (second  cycle)  was 
calculated  as  the  percentage  ratio  of  the  final  to  initial 
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dry  sample  weights. 


SWELLING  STRAIN  INDEX 

Unconfined  swelling  tests  were  carried  out  on  cylindrical 
samples  of  1.6  and  2.6  inch  diameter.  The  samples  were 
placed  in  a  metal  container  with  the  flat  faces  horizontal 
A  dial  gauge  was  positioned  on  the  top  surface  to  measure 
the  vertical  expansion.  The  specimen  was  submerged  in  tap 
water  at  room- temperature,  and  readings  were  taken  at 
i nterva 1 s . 

The  swelling  strain  index  (unconfined)  was  calculated 
from  the  following  equation  (Gyenge  <£  Herget,  1977): 


SSI  =  100. d/L 


SSI  =  Swelling  Strain  Index  (unconfined  %) 
d  =  Maximum  vertical  expansion  (mm) 

L  =  Specimen  vertical  dimension  (mm) 
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SILKSTONE  GEOTECHNICAL  TESTING  -  TRIAXIAL 


SAMPLE 

TYPE 

MOISTURE  % 

(psi) 

°3 

(psi) 

REMARKS 

J5A 

SANDSTONE 

4.0 

10010 

500 

SHEARS  2  5°  TCA 

J5B 

SANDSTONE 

12700 

1000 

- 

J4A 

SANDSTONE 

- 

6150 

0 

- 

J4B 

SANDSTONE 

- 

5600 

0 

Eo'  0 . 85xl0^psi 

J9A 

SANDSTONE 

2.6 

14950 

1500 

- 

J9B 

SANDSTONE 

2.8 

12470 

1000 

SHEAR  ~  2  8°  TCA 

J10A 

SANDSTONE 

- 

16630 

2000 

- 

J10B 

SANDSTONE 

3.5 

17930 

3000 

SHEAR  -  33°  TCA 

J14  A 

MUDSTONE 

5.0 

7820 

0 

COMPLICATED 
CONE  FAILURE 

J14B 

MUDSTONE 

00 

• 

9420 

500 

SHEAR ~  2  5°  TCA 

J16A 

MUDSTONE 

4.2 

13040 

1500 

SHEAR  ~  20°  TCA 

J17A 

MUDSTONE 

4.0 

15670 

2000 

COMPLICATED 

J17B 

MUDSTONE 

4 . 0 

11790 

1000 

SHEAR  ~  20°  TCA 

C10A 

SILTSTONE 

3.1 

12080 

1000 

SHEAR ~  2  8°  TCA 

C10B 

SILTSTONE 

3.4 

12040 

500 

CONE  FAILURE 

F7 

COAL 

- 

10140 

1000 

SHEAR-'  2  0°  TCA 

G14A 

COAL 

- 

5350 

500 

SHEAR  — '  2  8°  TCA 

G14B 

COAL 

- 

2270 

0 

COMPLICATED 

G1 3  A 

COAL 

- 

15690 

2500 

SHEAR ~  2  5°  TCA 

F6B 

COAL 

- 

12820 

1500 

COMPLICATED 

F3 

COAL 

- 

14710 

2000 

SHEAR ~ 2 5° TCA 

*  Samples  undercored  to  1.6  inch  diameter  along  core  axis. 
All  bedding  dips  between  10°  and  15°  TCA  (To  Core  Axis) 
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SILKSTONE  GEOTECHNICAL  TESTING  -  BRAZILIAN  DISC 


SAMPLE 

TYPE 

TENSILE 

STRENGTH 

(psi) 

REMARKS 

J4C 

Sandstone 

295 

- 

J9C 

Sandstone 

308 

- 

J9E 

Sandstone 

269 

- 

J9G 

Sandstone 

290 

— 

J5D 

Sandstone 

335 

- 

J10C 

Sandstone 

357 

- 

J10E 

Sandstone 

355 

— 

C9 

Coal 

119 

- 

G14C 

Coal 

97 

— 

Cl  IB 

Si ltstone 

588 

- 

C14B 

Siltstone 

389 

- 

J18B 

Siltstone 

626 

— 

C12A 

Mudstone 

864 

Possible  bad  failure 

F1C 

Mudstone 

282 

Weakness  plane 

J17D 

Mudstone 

197 

- 

J18D 

Mudstone 

631 

- 

X2A 

Mudstone 

1333 

Multiple  failure 

C15B 

Mudstone 

891 

- 

J4D 

Sandstone 

434 

Moisture  3.0% 

J9D 

Sandstone 

305 

Moisture  3.3% 

J9F 

Sandstone 

262 

- 

J5C 

Sandstone 

471 

Moisture  4.7% 

J5E 

Sandstone 

309 

- 

J10D 

Sandstone 

289 

- 

J10F 

Sandstone 

362 

Moisture  4.1% 

F6A 

Coal 

48 

- 

C13 

Coal 

263 

- 

C14A 

Siltstone 

304 

- 

G19B 

Siltstone 

364 

- 

J18C 

Siltstone 

745 

- 

FIB 

Mudstone 

375 

- 

J17C 

Mudstone 

457 

- 

J17E 

Mudstone 

171 

- 

XlB 

Mudstone 

96 

Complex  failure 

X2B 

Mudstone 

1043 

— 

*  Axis  of  discs  all  parallel  to  core  axis. 

*  Determination  of  tensile  test  plane  in  relation 
to  bedding  was  not  possible. 
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SILKSTONE  GEOTECHNICAL  TESTING  -  DIRECT  SHEAR 


SAMPLE 

TYPE 

AREA 

cm1 

NORMAL 

LOAD 

(N) 

FRICTION 

ANGLE 

(DEG) 

DESCRIPTION 

J3A 

SANDSTONE 

27 

950 

30.2 

SAW  CUT  PLANE  AT 

4  0°  TCA. 

ROCK  POWDER  LEFT 
AFTER  SHEARING 

1550 

31.8 

2550 

32.5 

3450 

33.1 

J12 

MUDSTONE 

30 

1350 

20.9 

JOINT  ~4  0°  TCA 

MEDIUM  ROUGH  WITH 
SOME  WAVINESS  ALONG 
STRIKE.  FLAKE 

DEBRIS  AFTER  SHEAR 

2850 

20.0 

4350 

18 . 7 

G10 

MUDSTONE 

38 

1350 

26.5 

JOINT  ^6  0°  TCA 
SLIGHTLY  ROUGH, WAVY 
ALONG  STRIKE  &  DIP. 
MAJOR  PORTION  FAILED 
AT  2850N . 

2850 

27.7 

G5 

MUDSTONE 

31 

1350 

21.7 

BEDDING  ~  75°  TCA 

V . THIN  CARBONACEOUS 
LENSES.  ROUGH  AND 
WAVY.  MUCH  SHEAR 
DEBRIS  AFTER  TEST. 

2850 

24 . 8 

4350 

23.3 

C2 

COAL 

36 

1350 

23.2 

BEDDING  ~  20°  TCA 

WAVY  ALONG  STRIKE, 
ROUGH  WITH  SOME 
SMOOTH  PATCHES. 

MUCH  SHEAR  DEBRIS 
AFTER  TEST. 

2850 

23.2 

4350 

25.1 

*  Residual  values  are  presented ,  as  there  were  no 
significant  peak  values. 

*  All  shearing  was  along  the  discontinuity  dip. 
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SILKSTONE  GEOTECHNICAL  TESTING  -  SLAKE  DURABILITY  INDEX 


SAMPLE 

TYPE 

INDEX  % 

REMARKS 

C7 

MUDSTONE 

78.5 

MUCH  FINE  BLOCKY  SEDIMENT 

G18 

MUDSTONE 

16.8 

MUCH  FINE  BLOCKY  SEDIMENT 

J15 

MUDSTONE 

95.3 

SOME  SANDY  SEDIMENT 

*  Second  cycle  index  calculated. 

*  Tap  water  at  20°  C  used. 


SILKSTONE  GEOTECHNICAL  TESTING  -  UNCONFINED  SWELLING 


SAMPLE 

TYPE 

DURATION 

HOURS 

SSI  % 

REMARKS 

C15B 

SILTSTONE 

72 

2.9 

SAMPLE  BROKEN  ALONG  BEDS 
WITH  CARBONACEOUS  LENSES 

J18A 

MUDSTONE 

18 

2.5 

COMPLETE  DISINTEGRATION 
OF  SAMPLE  AFTER  2  DAYS 

*  Tap  water  at  20° C  used. 

=  Swelling  Strain  Index  (unconfined) 


*  SSI 


Axial  Stress  (p si) 


Sandstone  G7 


IX  Ooo- 


II  OOo  - 


lo  OOO  - 


9  Ooo  " 


8  Ooo  ~ 


7ooo- 


CxPLOSive  g*iTTi_E 
PPILOO.E 


5  Ooo 


So o0  - 


•hOOO- 


I  ooo  - 


2ooo  - 


I  ooo -  1 


looo 


£  *10 


M I C  ROSTRA  I  N 


Strain  g-aug-ed  uniaxial  compRessioN 


Axial  Stress  cpsi) 
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Sandstone  T7A 


i 


lo  ooo  - 


A  ooo  - 


8  ooo  - 


7  OOO  - 


6  ooo- 


5  ooo  - 


4-  ooo  - 


3  ooo  - 


2.  OOo  - 


I  OOO- 


/ 


/ 


/• 

’  \ 

BRITTLE  FAILURE 


/vr 


L  ATERRL 
CMe<^A*riVE^ 


/ 


X 

/ 

X 


/ 

X 


Hd  -  512.*  IO* 


flVIflL 

E  =  i-u-*  lo1 


X 

I 

X 

I 

X 

I 

X 

I 

X 


/ 


/ 


I 

X 

I 

lx 

I 

X 


/ 


/<  =  Poisscn's  Uprno  «= —  =  0-27  4- 


I 


/ 


/ 


/ 


/ 


■ 


/ 


/ 


— nr - - — i - 1 - — i  i 

2000  LOO o  6ooo  goo©  loooo 

niCRO  STRAIN 


£«  io< 


Strain  Gauged  Uniaxial  ComPRession 


Axial  Stress  (p si) 
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'-'fi 
II  OOo  - 

lo  OOO  - 

^  OOO  - 

8  ooo  - 

7  OOO  - 

6  ooo 

Sooo  - 

4-  ooo  - 

3  ooo  - 


2 ooo  -| 


I  ooo- 


Sandstone  T  7B 


/ 


/i 


J  BRITTLE  FAH-OXE 


x 

/ 

X 


/ 

x 


/ 

X 

/ 


LATERAL 
(  NE&ftTive  ) 

X  IO* 


flxiAL 

E  =  Z*  I  6  X  i  o£ 


/ 


/ 


x 

I 

X 


/ 


X 

/ 

X 

I 

X 

I 

X 

I 

X 

I 

X 

I 

X 


/ 


/ 


/ 


/ 


/ 


/ 


/ 


X 

J 

X 


/ 


LX  =  Poisson's  &RTIC  =  JL  =  031 


/ 


/ 


./ 


/ 


T 


T 


T 


T 


T 


X 


I  OOO  2.000 


sooo  4ooo  Sooo  6000 
n  I  C  ROSTRA | N 


7000  e*io6 


Strain  G-auG-Erp  Uniaxial  Compression 


Axial  Stress  Cpsz's 
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Siltstone  Cl5R 


Strain  Gauged  Uniaxial  Compression 


■ 
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O 


UJ 

2 

o 

h 

CO 

CO 


\ 


ttJ 

K 

a 

a 

c 
a 

M 

J 

« 

%  « 

\ 

*\ 


\ 


V 


! 


-j 

<r 

x 

<x 


vO 

O 

X 

un 

00 

ii 

LU 


to 

00 

6 

II 

ui|  2 

II 

o 


z 

c> 

V/i 

w 

<£ 


it 

CC  V- 

IM  <t 

w  * 

51 

vO 


X 

O 

II 

u 

Ui 


\ 


\D 

O 

j 


\ 


'*■«. 


K*  In. 


V 

\ 


— I - - 


0 

0 

Q 

O 


O 

o 

o 

J- 
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o 

o 

o 

d 


“I 

o 

o 

o 

o 


o 

o 

o 

GO 


~ T~ 

O 

o 

o 

u> 


~~r 

o 

o 

o 

> 


o 

o 

o 

d 


cis*;  ss^^j-s  lu,xtf 


loeo  iooo  sooo  4-ooo  5ooo  6oo©  7  ooo  Sooo  l»oo 

pIlCROSTRfllN 

Strain  GfluirEP  Uniaxial  C  ornpR  eSSi  ON 


A* i«l  Stress  (pst  *  io~3) 
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II- 


10- 


S- 


8- 


7- 


6- 


if- 


3- 


Mudstone  F5 


x 

/ 

X 


X 

/ 


•'  \ 

f 


BRITTLE 
PAIL  U*E 


LflTETRR  L. 


T/  V£  ) 


=  7  Z8  x  lo 


/ 

jr 

f 

/ 

jc 


Y 

X 

I 

V 
0 

r 

X 


I 

X 

I 

X 

I 

X 

# 

IT 


A  XlftL 
E  =  115 


If 

I 

X 

I 

If 

» 

X 

I 

« 

I 

* 

# 

X 

t 

X 

I 

X 


/ 


^  ~  Poisson’s  Krtic  =  J=-  =  0-172 


1-i 


/ 


y 


I  ooo 


zooo 


ZOO  O 


if  ooo 


s  ooo 


6  ooo 


7  OOO 


8  ooo 


mi CROSTftAI N 


e*  io 


- r* 

9  OOO 

6 


Strain  Gt/aug-ed  (Jniaxial  Cqpipre  ss ion 


APPENDIX  8 


MINE  SIMULATION  RESULTS 
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»W»  IISH.  <,  »i$n..  USCOVTINVITT 


t 

i  S  5  •  • - 5 - i 

mirmfL  urrwer 

?.  mm  iisn.  «  »isn.  DiSBcnniKum 


jawiniuisMt 


RUN 

21  COAL 

VALLET 

CONFIGU 

RATION 

N«m= 

1 

En=* 

39E+04 

Vn=.19 

80 

Ett=- 

39E+04 

V  n=  •  1  9 

HDmc 

=  0.50 

En=- 

39E+04 

V  n= •  19 

WT»m 

=  0.025 

&n=- 

1 6E+0U 

S1G„ 

/S1Gtt= 

1 .000 

SM« 

^  KM 

6mm 

THICK 

DEPTH 

1  . 390E+03  . 1 50E*03  2.50  150.0 


ONSEfiM:  seam  1  segs  i-80 
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280 


.  ■  t- 

.  •  tw 

*  ■ 


RUN  21  COAL 

VALLEY  CONFIGURATION  2 

N  uik:  1 

En=.  39E*-04  Vn*.  19 

Nm.=  80 

E  TT= .  3  9  E  ♦  0  4  VfJ=  .  19 

WD„LpC  0.50 

Elt=  .  39E+04  V n=  .  19 

WW  0.025 

Gn=. 16E+04 

SIGn/SIGtr= 

1 . 000 

SH«s 

G«*.  thick  oepth 

1  .390E+03 

.  1 50E  +  03  2.50  150.0 

OFFSERM: 

Y  =  -14  5.0  SEGS  1-80 

*  •  »I6. 
.  •  si'. 
.  ■  sit. 


SHEAR  <  HORHflL  OISPlflCEHENTS 


— i - r - » — 

«  «  hi  i  «  a 

_ miyrflLjimg.. 
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T 


RUN  21  CORL  VRLLET  CONFIGURATION  2 


E„=. 39E+04 

V„=.19 

N  in  =  80 

E„=. 39E+0U 

CD 

• 

II 

c 

;> 

WD^,=  0.50 

E„=. 39E*0U 

Vn=.19 

NTMn=  0.025 

Gn=. 1 6E+01 

SIG«/SlGn-  1 

.  000 

SM* 

GiUM  THICK 

DEPTH 

1  . 390E+03  . 

1  50E  +  03  2.50 

150.0 

S  3T 


"5  I  5  S  ~S  5  1 

wninm  Biyrwct _ _ 


1  3  I  5  3 

_ (anwnw,  tisi«p 


• - 3 - 3 


•ti* 


282 


t 

i 


283 


284 


285 


RUN  22  COAL 

VRLLET  CONFIGURATION  3 

W  1 

E„=. 39E+04 

Vw-.19 

Nftt=  60 

E  TT=  .  3  9  E  +  0  4 

V„=.19 

w°MLr=  i-oo 

E„- . 39E+04 

V  Tl=  -  1  9 

WTum=  0.025 

G„=. 1 6E  +  04 

SIG„/SIG„= 

1.000 

SM* 

THICK 

DEPTH 

1  .  3  9  0  E  0  3 

.  1  5  0  E  +  0  3  2.50 

150.0 

0  F  F  S  E  R  M  : 

Y  =  -145.0  SEGS  1-80 

J  T  t  A  1  N  J 


t. 


.  S16, 
'  *16. 

*  *16. 


SH  E.  A  R  »  NORMAL  DISPLACEMENTS 


HMi28inaL  hsirice 
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287 


288 


shear  displ.  «  01  spl.  uscDKimum 


t 

.• 


■wuwim.  Burner 

HO  A  HAL  D1 SPL.  *  D1SPL.  0 1 SC0H1 1  HU  1 T T 


289 


290 


RUN 

23  COAL 

VALLEY  CONFIGURATION  4 

N«m= 

1 

E  n-  .  3  9  E  +  0  4  V  1T=  .  1  9 

N«.= 

60 

ETT=  .  3 9 E  +  0  4  Vu=.19 

WD-l, 

=  0.50 

E  M=  .  3  9  E  +  0  4  V  Tl=  .  19 

WT„.n 

=  0.025 

GIt=.  16E+04 

S1G„ 

/sigtt= 

1 . 000 

S  M  * 

GKM  thick  depth 

1  . 

3  9  0  E  +  0  3 

. 1 5  0  E  +  0  3  2.50  1  50.0 

OFFSERM : 

Y  =  -145.0  SEGS  1-60 

SHEAR  4  NORMAL  DISPLACEMENTS 


.  *  «, 
.  •  U. 


P - 1 - ' 

-•  -m  -n 


ttMlIftflRL  DISIMg 
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S  T  I  ft  1  U  S 


i. 


sit. 

«»«. 

Sit. 


BUN  23  COAL 

VALLEY  CONFIGURATION  4 

N«."  1 

E  n=  .  3  9  E  +  0  4  V  lT=  .  19 

N„.=  80 

E„=.39E  +  04  V  M= .  19 

NDmlp=  0.50 

EIt=  .  3  9E  +  04  V TI=  .  1  9 

WTu>m=  0.025 

Git=.  16E  +  04 

sig„/sigtt= 

1  .000 

SM*  E.„. 

gmm  thick  depth 

1  .  3  9  0  E  +  0  3 

. 1  5  0  E  +  0  3  2.50  1  50.0 

OFFSEflM : 

Y  =  -148.0  SEGS  1-80 

SH£RI  t  KORHfil  01 SPLSCEHEMTS 


5  5  r  »  ■ 

itwiisnu  oisiMg 


s — s — 5  r  .  . 

_ MMiwnw,  umpire 
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